DQSpillFlags_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
///
/// \brief   Calculate the spill level data quality flags and 
///          save to sumdata::EventQuality
/// \authors xbbu@fnal.gov
/// \date    11/10/2014
///
////////////////////////////////////////////////////////////////////////

//C++ includes 
#include <cmath>
#include <map>
#include <string>
#include <vector>

// Framework 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_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"

// NOvASoft includes
#include "Calibrator/Calibrator.h"
#include "CMap/service/CMapService.h"
#include "DAQChannelMap/DAQChannelMap.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "SummaryData/EventQuality.h"
#include "RawData/DAQHeader.h"

namespace dqsf
{
  struct DCMId
  {
    DCMId(geo::View_t v, int _x, int _y) : view(v), x(_x), y(_y) {}
    geo::View_t view;
    int x, y;

    bool operator<(const DCMId rhs) const
    {
      return std::make_tuple(view, x, y) < std::make_tuple(rhs.view, rhs.x, rhs.y);
    }
  };

  class DQSpillFlags : public art::EDProducer 
  {
  public:
    explicit DQSpillFlags(fhicl::ParameterSet const& pset);
    virtual ~DQSpillFlags();

    void produce(art::Event& evt);

  protected:
    void FillCountVars(sumdata::EventQuality* spillqual,
                       const std::vector<rb::CellHit>& chits) const;
    void FillSyncVars(sumdata::EventQuality* spillqual,
                      const std::vector<rb::Cluster>& slices,
                      const std::vector<rb::Track>& tracks) const;
    void FillRawVars(sumdata::EventQuality* spillqual,
                     const rawdata::DAQHeader& daqheader);
    
    void FillPopulatedDCMs(const rb::Cluster& clust,
                           std::set<DCMId>& pop) const;
    
    void AccumulateSyncMetric(const rb::Cluster& clust,
                              const std::set<DCMId>& pop,
                              int& nedge, int& nmatch) const;
    
    std::string  fCellHitLabel; ///< Label for cell hits
    std::string  fSliceLabel;   ///< Label for slices
    std::string  fTrackLabel;   ///< Label for tracks
    std::string  fRawDataLabel; ///< Label for raw data
    
    bool fTrackSyncMetric;      ///< Use track hits rather than all hits
    
    double       fTimeLow;      ///< Hit time of beam trigger window low end
    double       fTimeHigh;     ///< Hit time of beam trigger window high end
    unsigned int fNoisyHitMin;  ///< Miniumum number of hits to be considered a noisy APD
  };
  
  
  //..............................................................
  DQSpillFlags::DQSpillFlags(fhicl::ParameterSet const& pset)
  : EDProducer(pset)
  {
    fCellHitLabel = pset.get< std::string >("CellHitLabel");
    fSliceLabel   = pset.get< std::string >("SliceLabel");
    fTrackLabel   = pset.get< std::string >("TrackLabel");
    fRawDataLabel = pset.get< std::string >("RawDataLabel");
    fTrackSyncMetric = pset.get< bool     >("TrackSyncMetric");
    fTimeLow      = pset.get< double      >("TimeLow");
    fTimeHigh     = pset.get< double      >("TimeHigh");
    fNoisyHitMin  = pset.get< int         >("NoisyHitMin");
    
    produces< sumdata::EventQuality >();
  }
  
  //......................................................................
  DQSpillFlags::~DQSpillFlags()
  {
  }
  
  //......................................................................
  void DQSpillFlags::produce(art::Event& evt)
  {
    std::unique_ptr<sumdata::EventQuality> spillqual(new sumdata::EventQuality);

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

    //get all hits info         
    art::Handle<std::vector<rb::CellHit> > chits;
    evt.getByLabel(fCellHitLabel, chits);

    art::Handle<std::vector<rb::Cluster>> slices;
    art::Handle<std::vector<rb::Track>> tracks;

    art::Handle<rawdata::DAQHeader> daqheader;
    evt.getByLabel(fRawDataLabel,daqheader);
    if( !daqheader.failedToGet() )
      FillRawVars(spillqual.get(),*daqheader);

    FillCountVars(spillqual.get(), *chits);

    //skip caclulating FD track sync variables
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      evt.put(std::move(spillqual));
      return; //skip caclulating FD track sync variables
    }
 
    evt.getByLabel(fSliceLabel, slices);

    if(fTrackSyncMetric) evt.getByLabel(fTrackLabel, tracks);

    if(fTrackSyncMetric)
      FillSyncVars(spillqual.get(), *slices, *tracks);
    else
      FillSyncVars(spillqual.get(), *slices, {});
                        
    
    //save info to event tree
    evt.put(std::move(spillqual));
    
  }//end event loop

  //..............................................................
  void DQSpillFlags::FillCountVars(sumdata::EventQuality* spillqual,
                                   const std::vector<rb::CellHit>& chits) const
  {
    const daqchannelmap::DAQChannelMap* cmap = art::ServiceHandle<cmap::CMap>()->Map();
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<geo::LiveGeometry> livegeom;

    unsigned int horizhits = 0; // All hits for horizontal modules
    unsigned int ndcm3hits = 0; // DCM3 hits for horizontal modules

    unsigned int ncells_outtime = 0; // # of cells in out-beam-windomw

    // Chosen an arbitrary diblock/dcm/feb/pixel to encode a dchan in
    // order to get diblock type (instead of hardcoding). This allows
    // one to set the max total dcm value for either ND or FD
    const daqchannelmap::dchan dchan = cmap->encodeDChan(geom->DetId(), 1, 1, 0, 0);
    const daqchannelmap::DiBlock_TYPE db_t = cmap->getDiBlockType(dchan);
                                                        
    // Irrespective of ND Muon Catcher, this sets max number of DCM within 
    // detector diblock to set limits on loop over DCMs later on
    const unsigned int totdcm = cmap->getNumberOfDCM(db_t);
    const unsigned int totdb  = cmap->getTotalNumberOfDiblocks();

    std::vector<std::vector<unsigned int> > ndcmhits(totdb);
    for(auto & itr : ndcmhits) itr.resize(totdcm, 0);

    std::map< unsigned int, 
              std::map<unsigned int, unsigned int> > dbapdmap; // Map of diblock/APD/nhits

    for(const rb::CellHit& chit: chits){

      const unsigned int dcm = cmap->getDCM(chit.DaqChannel());
      const unsigned int db  = cmap->getDiBlock(chit.DaqChannel());
      const unsigned int apd = cmap->getFEB(chit.DaqChannel());
      // Offset by 1 as db and dcm accounting starts at 1
      ++ndcmhits[db-1][dcm-1];

      // Only look at hits outside of beam window
      if(chit.TNS() < fTimeLow || chit.TNS() > fTimeHigh ){
        ++ncells_outtime; // total number of cells out of time
        
        // Cells 95, 87, 79 and 71 in Y-view in ND are impacted by
        // lights on the catwalk. See K. Sachdev doc-db 12426
        if(geom->DetId() == novadaq::cnv::kNEARDET &&
           chit.Plane() < geom->FirstPlaneInMuonCatcher() &&
           chit.View() == geo::kY){
          ++horizhits;
          int cell = chit.Cell();
          if(cell == 95 || cell == 87 || cell == 79 || cell == 71)
            ++ndcm3hits;
        } // DB horizontal plane only, excluding Muon Catcher
        
        // count number of hits for each apd into ndbapd map
        ++dbapdmap[db][apd];
      } // outside beam window
    } // end loop over all hits
    
    unsigned int ndeaddcms = 0;
    
    // Loop over di-blocks
    for(unsigned int dbIdx = 0; dbIdx < totdb; ++dbIdx){
      if(geom->DetId() == novadaq::cnv::kNEARDET && (dbIdx < (totdb-1))){
        for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
          if(ndcmhits[dbIdx][dcmIdx] == 0) ndeaddcms++;
        }
      } // end loop for ND minus  Muon Catcher
      else if(geom->DetId() != novadaq::cnv::kNEARDET){
        for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
          if(ndcmhits[dbIdx][dcmIdx] == 0) ++ndeaddcms;
        }
      } // end loop for NOT ND
    }  // end loop over diblocks
    // Special case for the 2 Near Det Muon Catcher DCMS
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      for(unsigned int dcmIdx = 0; dcmIdx < 2; ++dcmIdx){
        if(ndcmhits[totdb-1][dcmIdx] == 0) ++ndeaddcms;
      } // end loop for ND Muon Catcher DCMs
    } // end loop for ND
    
    // Calculate fraction of partial DCM3 hits in horizontal modules to
    // detect lights-on effect in ND
    float filter_lightson = 0.; // fraction of partial DCM3 hits in horizontal modules
    if(horizhits > 0) filter_lightson = float(ndcm3hits)/horizhits;
    
    // Noisy APD if it has at least fNoisyHitMin hits.
    // Original value of fNoisyHitMin = 28 for ND
    // See X. Bu doc-db 12316
    // Total number of noisy apd in ND
    int napd_hitmin_tot = 0;
    
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      for(const auto dbItr : dbapdmap){
        for(const auto apdItr : dbItr.second){
          if(apdItr.second > fNoisyHitMin) napd_hitmin_tot++;
        } // end loop APD/Value pair
      } // end loop dbapdmap
    } //end detector check
    
    //save the values for dq flags
    spillqual->nmissingdcmslg = livegeom->NDropouts();
    spillqual->nmissingdcms   = ndeaddcms;
    spillqual->fracdcm3hits   = filter_lightson;
    spillqual->nouttimehits   = ncells_outtime;
    spillqual->nnoisyapds     = napd_hitmin_tot;
  }

  //..............................................................
  void DQSpillFlags::FillSyncVars(sumdata::EventQuality* spillqual,
                                  const std::vector<rb::Cluster>& slices,
                                  const std::vector<rb::Track>& tracks) const
  {
    // Which DCMs are populated at all, so we don't penalize at edges and for
    // missing DCMs (that's a different cut)
    std::set<DCMId> pop;

    for(const rb::Cluster& slice: slices) FillPopulatedDCMs(slice, pop);

    // How many edge crossings we considered
    int nedge = 0;
    // How many had a match
    int nmatch = 0;

    if(fTrackSyncMetric){
      for(const rb::Track& track: tracks)
        AccumulateSyncMetric(track, pop, nedge, nmatch);
    }
    else{
      for(const rb::Cluster& slice: slices){
        if(!slice.IsNoise()) AccumulateSyncMetric(slice, pop, nedge, nmatch);
      }
    }
    
    spillqual->dcmedgematchfrac = (nedge ? nmatch/double(nedge) : -1);
  }

  enum Dir{kLeft, kUp, kRight, kDown};

  //..............................................................
  Dir Opposite(Dir d)
  {
    switch(d){
    case kLeft:  return kRight;
    case kRight: return kLeft;
    case kUp:    return kDown;
    case kDown:  return kUp;
    default: assert(0 && "Not reached");
    }
  }

  //..............................................................
  DCMId Move(DCMId dcm, Dir d)
  {
    switch(d){
    case kLeft:  return DCMId(dcm.view, dcm.x-1, dcm.y  );
    case kRight: return DCMId(dcm.view, dcm.x+1, dcm.y  );
    case kUp:    return DCMId(dcm.view, dcm.x,   dcm.y+1);
    case kDown:  return DCMId(dcm.view, dcm.x,   dcm.y-1);
    default: assert(0 && "Not reached");
    }
  }

  struct Edges
  {
    Edges() : flag() {}
    bool flag[4];
  };

  //..............................................................
  void DQSpillFlags::FillPopulatedDCMs(const rb::Cluster& clust,
                                       std::set<DCMId>& pop) const
  {
    for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
      art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);

      const int dcmx = chit->Plane() / 64;
      const int dcmy = chit->Cell() / 64;

      pop.insert(DCMId(chit->View(), dcmx, dcmy));
    }
  }

  //..............................................................
  void DQSpillFlags::AccumulateSyncMetric(const rb::Cluster& clust,
                                          const std::set<DCMId>& pop,
                                          int& nedge, int& nmatch) const
  {
    std::map<DCMId, Edges> edges;

    for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
      art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);

      const int modplane = chit->Plane() % 64;
      const int modcell = chit->Cell() % 64;

      const int dcmx = chit->Plane() / 64;
      const int dcmy = chit->Cell() / 64;

      const DCMId key(chit->View(), dcmx, dcmy);

      if(modplane == 0 || modplane == 1) edges[key].flag[kLeft] = true;
      if(modplane == 62 || modplane == 63) edges[key].flag[kRight] = true;
      if(modcell == 0) edges[key].flag[kDown] = true;
      if(modcell == 63) edges[key].flag[kUp] = true;
    } // end for chitIdx

    // Make a copy because [] access to edges in the loop below might mutate
    // what we're looping over?
    std::map<DCMId, Edges> edges_it = edges;

    for(auto it: edges_it){
      const DCMId key = it.first;
      const Edges val = it.second;
      
      for(Dir dir: {kLeft, kRight, kUp, kDown}){
        if(!val.flag[dir]) continue;
        
        const DCMId nei = Move(key, dir);
          // Only count if the neighbouring DCM has hits at all
        if(!pop.count(nei)) continue;
        
        if(edges[nei].flag[Opposite(dir)]){
          ++nedge;
          ++nmatch;
        }
        else{
            // If there's no match, count two chances, because a good match
            // will have been counted from both sides.
          nedge += 2;
        }
      } // end for dir
    } // end for it
  }
  
  //-----------------------------------------------------
  void DQSpillFlags::FillRawVars(sumdata::EventQuality* spillqual,
                     const rawdata::DAQHeader& daqheader)
  {
     spillqual->nmicroslices = daqheader.TotalMicroSlices();
  }//end for FillRawVars
  
} // end namespace dqsf
////////////////////////////////////////////////////////////////////////

DEFINE_ART_MODULE(dqsf::DQSpillFlags)