AttenCache.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   Helper for AttenCurve
/// \author  Chris Backhouse - bckhouse@caltech.edu
///   
////////////////////////////////////////////////////////////////////////
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Calibrator/func/AttenCache.h"
#include "Calibrator/func/AttenCurve.h"
#include "CalibrationUtils/CalibUtil.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/GeometryBase.h"
#ifdef NOVACMAKE
#include "Database/cxx/include/Table.h"
#else
#include "Database/Table.h"
#endif

namespace calib
{
  //......................................................................
  AttenCache::AttenCache( bool is_data, int brightness_bin )
    : fUseEpochs(false),
      fReadEpochsFromCSV(false),
      fInitializedData(false), fInitializedMC(false),
      fInitializedUncalibratedMask(false),
      fGainSetting(0),
      fBrightnessBin(brightness_bin),
      fIsData(is_data)
  {
  }

  //......................................................................
  AttenCache::~AttenCache()
  {
    for(auto it: fMapMC)   delete it.second;
    for(auto it: fMapData) delete it.second;
  }

  //......................................................................
  void AttenCache::Add(AttenCurve* res)
  {
    if(fIsData){
      fMapData[Key_t(res->det, res->chan)] = res;
    }
    else{
      fMapMC[Key_t(res->det, res->chan)] = res;
    }
  }

  //......................................................................
  const AttenCurve* AttenCache::ForData(const geo::GeometryBase* geom,
                                        geo::OfflineChan chan)
  {
    // TODO - this won't behave right if we ask for more than one detector
    if(!fInitializedData){
      fInitializedData = true;
      LoadFromDatabase(geom);
    }

    return Get(TranslateChannelData(chan, geom));
  }

  //......................................................................
  const AttenCurve* AttenCache::ForMC(const geo::GeometryBase* geom,
                                      geo::OfflineChan chan, bool maskUncalibratedChannels)
  {
    // TODO - this won't behave right if we ask for more than one detector
    if(!fInitializedMC){
      fInitializedMC = true;
      LoadFromDatabase(geom);
    }

    if(!fInitializedUncalibratedMask && maskUncalibratedChannels){
      fInitializedUncalibratedMask = true;
      LoadFromDatabaseUncalibratedMask(geom);
    }

    if(maskUncalibratedChannels &&
       fMapUncalibratedMask.count(TranslateChannelData(chan, geom))){
      return 0;
    } 

    return Get(TranslateChannelMC(chan, geom));
  }

  //......................................................................
  const AttenCurve* AttenCache::Get(Key_t key) const
  {
    if(fIsData){
      auto it = fMapData.find(key);
      if(it != fMapData.end()) return it->second;
    }
    else{
      auto it = fMapMC.find(key);
      if(it != fMapMC.end()) return it->second;
    }

    return 0;
  }

  //......................................................................
  void AttenCache::GetEpochFromDB(int idet, bool forceData)
  {
    nova::dbi::Table tbl("Calibration/tables/CalibAttenEpochs.xml");
    if (fIsData || forceData)
      tbl.SetDataTypeMask(nova::dbi::kDataOnly);
    else
      tbl.SetDataTypeMask(nova::dbi::kMCOnly);
    tbl.SetTag(fEpochTag);
    tbl.SetDetector(idet);
    tbl.SetMinTSVld(fRun);
    tbl.SetMaxTSVld(fRun);
    tbl.SetVerbosity(100);
    assert(tbl.Load());
    assert(tbl.NRow() == 1);

    int epochIdx = tbl.GetColIndex("epoch");
    tbl.GetRow(0)->Col(epochIdx).Get(fEpoch);
  }


  
  //......................................................................
  std::string AttenCache::GetCSVPath(novadaq::cnv::DetId det,
                                     bool                points,
                                     bool                forceData) 
  {
    // necessary components of path
    const char* ups_path = getenv("CALIBCSVS_CSV_PATH");
    assert(ups_path);
    const char* pointsStr = points ? "points" : "consts";

    const char* detStr = calib::getDetString(det);
       //(det == novadaq::cnv::kFARDET) ? "fd" : "nd";
    const char* mcStr = (fIsData || forceData) ? "data" : "mc";

    // extra specifiers depending on epoch, gain, and brightness settings
    std::string gainStr, epStr, fbStr;
    fbStr = (forceData || fIsData || fBrightnessBin<0) ? "" : ".fb"+std::to_string(fBrightnessBin);

    //RJN check that directory exists
    std::string dirString=TString::Format("%s/%s",ups_path, fTag.c_str()).Data();
    if(!doesDirectoryExist(dirString)) {
       LOG_ERROR("AttenCache") << dirString << " is not a directory is there a mismatch between the CALIBCSVS_CSV_PATH and the calibration tag in Calibration.fcl ?\nCALIBCSVS_CSV_PATH is" << ups_path << "\nTag is " << fTag << "\n" ;
       abort();
       //       assert(!"Can't find the right calibration directory");
    }


    if(fReadEpochsFromCSV) {
      //Need to read the epochs and then assign things correctly
      //Currently this won't work unles fUseCSVsFromUPS is set RJN should fix this soon.
      std::string prePattern=TString::Format("%s/%s/calib_atten_%s.%s.%s.%s%s", 
                                             ups_path, fTag.c_str(),                  // -Path
                                             pointsStr, detStr, mcStr, fTag.c_str(),  // -CSV types
                                             gainStr.c_str()).Data();
      std::string postPattern=TString::Format("%s.csv",fbStr.c_str() ).Data();

      LOG_DEBUG ("AttenCache")   << "In GetCSVPath prePattern "
                                 << prePattern.c_str()
                                 << "\t postPattern " << postPattern.c_str() << "\n";
      
      return getCSVFilenameByParsingDirectory(fRun,prePattern,postPattern);

      
    }
    else if (fUseEpochs) {
      GetEpochFromDB(det,forceData);
      assert(!fEpoch.empty());
      epStr="."+fEpoch;   gainStr="";
    }
    else {
      epStr="";
      if (!fGainSetting)  gainStr="";
      else                gainStr=TString::Format(".g%d",fGainSetting);
    }


    // if configured to provide own path
    if(!fUseCSVsFromUPS){   
      return TString::Format("%s/calib_atten_%s%s.csv", 
                             fCSVSourceDir.c_str(),
                             pointsStr,
                             fbStr.c_str() ).Data();

    }


    //  Return final path                                |----|  last 3 could be empty
    return TString::Format("%s/%s/calib_atten_%s.%s.%s.%s%s%s%s.csv", 
                           ups_path, fTag.c_str(),                  // -Path
                           pointsStr, detStr, mcStr, fTag.c_str(),  // -CSV types
                           gainStr.c_str(), epStr.c_str(),          // -Extra specifiers which 
                           fbStr.c_str() ).Data();                  //    could be empty    
  }

  //......................................................................
  void AttenCache::LoadFromDatabase(const geo::GeometryBase* geom)
  {
    // For the FD this will make a cache ~80MB. That's OK, but denser LOWESS
    // points would increase this rapidly.

    LOG_VERBATIM("AttenCache") << "Initializing attenuation cache for "
                               << novadaq::cnv::DetInfo::GetName(geom->DetId())
                               << "\t fReadEpochsFromCSV:\t" << fReadEpochsFromCSV << "\n";

    const std::string constsCSVSource = GetCSVPath(geom->DetId(), false);
    LOG_VERBATIM("AttenCache") << "Loading attenuation cache consts from "
                               << constsCSVSource;

    FILE* fconsts = fopen(constsCSVSource.c_str(), "r");
    if(!fconsts){
      std::cerr << "Couldn't open " << constsCSVSource << std::endl;
      if( constsCSVSource.find(".fb") != std::string::npos )
        std::cerr << "   Note: cannot model fiber brightness while using"
                  << " UPS CVS constants tagged v9 or before." << std::endl;
      abort();
    }

    char hdr[4096];
    fscanf(fconsts, "%s", hdr);
    // Expect the file to start with this exact header
    if(hdr != std::string("#coeff_exp,atten_length,background,edge_low,edge_high,coeff_low,coeff_high,chisq,channel,tv")){
      std::cerr << "Bad header in " << constsCSVSource
                << " : '" << hdr << "'" << std::endl;
      abort();
    }

    unsigned int line=0;
    while(!feof(fconsts)){
      line++;
      AttenCurve* res = AttenCurve::Uninitialized(geom->DetId(),
                                                  geo::OfflineChan(-1, -1));

      int channel, tv;
      const int n = fscanf(fconsts, "%f,%f,%f,%f,%f,%f,%f,%f,%d,%d",
                           &res->coeff_exp, &res->atten_length, &res->background,
                           &res->edge_low, &res->edge_high,
                           &res->coeff_low, &res->coeff_high,
                           &res->chisq,
                           &channel, &tv);
      if(n == EOF){
        delete res;
        break;
      }

      if(n != 10){
        std::cerr << "Bad row in " << constsCSVSource
                  << "\n  Read " << n << " columns on line " 
                  << line+1 << std::endl;
        abort();
      }

      res->chan = geo::OfflineChan::FromDBValidityChan(channel);

      bool inMuonCatcher = false;
      const geo::PlaneGeo* geoplane = 0;

      const int plane = res->chan.Plane();
      if(fIsData && plane < 1000){
        // Normal data just has regular channel numbers
        geoplane = geom->Plane(plane);
        if(geom->DetId() == novadaq::cnv::kNEARDET &&
           plane >= int(geom->FirstPlaneInMuonCatcher())){
          inMuonCatcher = true;
        }
      }
      else{
        // Otherwise is merged views or MC
        geoplane = RepresentativePlane(geom, plane, inMuonCatcher);
      }

      if(!geoplane) continue;
      const geo::CellGeo* geocell = geoplane->Cell(res->chan.Cell());
      if(!geocell) continue;

      res->cell_length = 2*geocell->HalfL();

      if(inMuonCatcher){
        // Find a plane from the same view that's not in the muon catcher, to
        // find out how much the muon catcher cell must be displaced.
        int i = 0;
        const geo::PlaneGeo* regularPlane = geom->Plane(i);
        while(regularPlane->View() != geoplane->View())
          regularPlane = geom->Plane(++i);
        const geo::CellGeo* regularCell = regularPlane->Cell(res->chan.Cell());

        // Assume that the cells are aligned on the readout side
        res->center_offset = regularCell->HalfL()-geocell->HalfL();
      }

      res->initialized = true;

      Add(res);

      // Figure out which diblock, mark it calibrated
      const int db = plane/64+1;
      fDiblocks.insert(db);
    } // end while

    fclose(fconsts);


    // Now fetch the LOWESS points and add them to the curves we've just
    // created.

    const std::string pointsCSVSource = GetCSVPath(geom->DetId(), true);
    mf::LogInfo("AttenCache") << "Loading attenuation cache points from "
                              << pointsCSVSource;

    FILE* fpoints = fopen(pointsCSVSource.c_str(), "r");
    if(!fpoints){
      std::cerr << "Couldn't open " << pointsCSVSource << std::endl;
      abort();
    }

    // Expect the file to start with this exact header
    fscanf(fpoints, "%s", hdr);

    // Expect the file to start with this exact header
    if(hdr != std::string("#w,factor,channel,tv")){
      std::cerr << "Bad header in " << pointsCSVSource
                << " : '" << hdr << "'" << std::endl;
      abort();
    }

    int totPoints = 0;
    line = 0;
    while(!feof(fpoints)){
      line++;
      float w, factor;
      int channel, tv;

      const int n = fscanf(fpoints, "%f,%f,%d,%d", &w, &factor, &channel, &tv);
      if(n == EOF) break;

      if(n != 4){
        std::cerr << "Bad row in " << pointsCSVSource
                  << "\n  Read " << n << " columns on line " 
                  << line+1 << std::endl;
        abort();
      }

      ++totPoints;

      const geo::OfflineChan chan = geo::OfflineChan::FromDBValidityChan(channel/1000);

      AttenCurve* curve;

      if(fIsData){
        auto it = fMapData.find(Key_t(geom->DetId(), chan));
        if(it == fMapData.end()){
          std::cerr << "Couldn't find appropriate curve for atten point. "
                    << "Chan = " << chan << std::endl;
          std::abort();
        }
        curve = it->second;
      }
      else{
        auto it = fMapMC.find(Key_t(geom->DetId(), chan));
        if(it == fMapMC.end()){
          std::cerr << "Couldn't find appropriate curve for attenp point. "
                    << "Chan = " << chan << std::endl;
          std::abort();
        }
        curve = it->second;
      }

      const unsigned int interpIdx = channel%1000;
      if(curve->interp_pts.size() <= interpIdx)
        curve->interp_pts.resize(interpIdx+1);
      curve->interp_pts[interpIdx] = {w, factor};
    } // end while

    fclose(fpoints);

    // Save a little memory
    for(auto it: fMapData) it.second->interp_pts.shrink_to_fit();
    for(auto it: fMapMC) it.second->interp_pts.shrink_to_fit();

    // Check that the number of points we got from the DB is the same as the
    // total in all the AttenCurves. If these numbers don't match, it means
    // there were duplicates or gaps.
    int pointsCount = 0;
    if(fIsData){
      for(auto it: fMapData) pointsCount += it.second->interp_pts.size();
    }
    else{
      for(auto it: fMapMC) pointsCount += it.second->interp_pts.size();
    }

    if(pointsCount != totPoints){
      std::cerr << "Mismatch between DB and AttenCurve counts. "
                << pointsCount << " vs "<< totPoints << std::endl;
      std::abort();
    }
  }

  //......................................................................
  void AttenCache::LoadFromDatabaseUncalibratedMask(const geo::GeometryBase* geom)
  {
    mf::LogInfo("AttenCache")
      << "Initializing data attentuation fit chisq values for " 
      << novadaq::cnv::DetInfo::GetName(geom->DetId()) 
      << " from database to mask uncalibrated channels...\n";
    
    nova::dbi::Table tblConsts("Calibration/tables/CalibAttenConstsVld.xml");
    tblConsts.SetDetector(geom->DetId());
    tblConsts.SetDataTypeMask(nova::dbi::kDataOnly);

    if(!fTag.empty()) tblConsts.SetTag(fTag);
    tblConsts.SetVerbosity(100);

    tblConsts.SetMinTSVld(1);
    tblConsts.SetMaxTSVld(100);

    bool ok = true;

    if(!fCSVSourceDir.empty()){
      mf::LogInfo("AttenCache") 
        << "Trying to get data constants from the MC CSV file, "
        << "no good way to handle this if you want to use masking. Aborting." 
        << std::endl;
      std::abort();
    }
    else{
      if(fUseCSVsFromUPS){
        // Always ask for data, because we're masking the MC to match real data
        // conditions.
        const std::string constsCSVSource = GetCSVPath(geom->DetId(), false, true);
        mf::LogInfo("AttenCache") << "Loading attenuation cache from "
                                  << constsCSVSource << std::endl;
        ok = ok && (tblConsts.LoadFromCSV(constsCSVSource) > 0);
        mf::LogInfo("AttenCache") << "Got " << tblConsts.NRow() << " rows from CSV" << std::endl;
      }
      else{
        ok = ok && tblConsts.Load();
      }
    }

    if(!ok){
      // No safe way to handle this
      mf::LogError("AttenCache") << "Unable to load attenuation tables. Aborting"
                                 << std::endl;
      // I've definitely seen mf lose races with shutdown
      std::cerr << "Unable to load attenuation tables. Aborting"
                << std::endl;
      std::abort();
    }

    const int chisqCol = tblConsts.GetColIndex("chisq");

    const int Nconsts = tblConsts.NRow();
    for(int nconsts = 0; nconsts < Nconsts; ++nconsts){
      nova::dbi::Row* row = tblConsts.GetRow(nconsts);

      double chisq;
      ok = ok && row->Col(chisqCol).Get(chisq);

      if(!ok) std::abort();

      const geo::OfflineChan chan = geo::OfflineChan::FromDBValidityChan(row->Channel());

      if(chisq <= 0 || chisq > 0.2){
        fMapUncalibratedMask.insert(Key_t(geom->DetId(), chan));
      } 

      // Figure out which diblock, mark it calibrated
      const int db = chan.Plane()/64+1;
      fDiblocks.insert(db);
    } // end for nconsts
  }
  
  //......................................................................
  int AttenCache::GetView(int plane, const geo::GeometryBase* geom) const
  {
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      const int firstPlane = geom->FirstPlaneInMuonCatcher();
      const int isMuonCatcher = (firstPlane <= plane);
      return geom->Plane(plane)->View()+(2*isMuonCatcher);
    }

    return geom->Plane(plane)->View();
  }

  //......................................................................
  const geo::PlaneGeo* AttenCache::
  RepresentativePlane(const geo::GeometryBase* geom,
                      int view,
                      bool& inMuonCatcher) const
  {
    // Unfortunately a bunch of different conventions competing here

    // View-merged MC, or average constants data
    if(view == 0 || view == 1 || view == 1000 || view == 1001){
      inMuonCatcher = false;
      return geom->Plane(*geom->GetPlanesByView(geo::View_t(view)).begin());
    }

    // View-merged MC, or some scheme I'm not sure is still in use, or
    // averaged-constants MC
    if(view == 2 || view == 3 || view == 10 || view == 11 ||
       view == 1002 || view == 1003){

      if(geom->DetId() != novadaq::cnv::kNEARDET){
        std::cerr << "Muon catcher constants in non-ND detector" << std::endl;
        std::abort();
      }

      inMuonCatcher = true;
      view = view%2;

      const int mcp = geom->FirstPlaneInMuonCatcher();
      const geo::PlaneGeo* geoplane = geom->Plane(mcp);
      if(geoplane->View() != geo::View_t(view))
        geoplane = geom->Plane(geom->NextPlaneOtherView(mcp, +1));

      return geoplane;
    }

    std::cerr << "Bad encoded view " << view << std::endl;
    std::abort();
  }

  //......................................................................
  AttenCache::Key_t AttenCache::
  TranslateChannelData(geo::OfflineChan chan,
                       const geo::GeometryBase* geom) const
  {
    const int db = chan.Plane()/64+1;
    if(!fDiblocks.count(db)){
      // No calibration for this diblock, use the averages instead.
      // Stored at channel 1,000,000+1000*view+cell
      const int view = GetView(chan.Plane(), geom);
      return Key_t(geom->DetId(),
                   geo::OfflineChan(1000+int(view), chan.Cell()));
    }

    return Key_t(geom->DetId(), chan);
  }

  //......................................................................
  AttenCache::Key_t AttenCache::
  TranslateChannelMC(geo::OfflineChan chan,
                     const geo::GeometryBase* geom) const
  {
    // If it's a muon catcher plane, check to see if there is a special entry
    // for it
    if(geom->DetId() == novadaq::cnv::kNEARDET &&
       chan.Plane() >= int(geom->FirstPlaneInMuonCatcher())){
      const int viewInt = GetView(chan.Plane(), geom);
      const Key_t ret(geom->DetId(), geo::OfflineChan(viewInt, chan.Cell()));
      if(Get(ret)) return ret;
    }
    // If not, fall back to the non-muon-catcher constants
      
    // This is the way we encode the MC parameters. Plane represents the view.
    // TODO, we should move to plane+1000 like the data?
    return Key_t(geom->DetId(),
                 geo::OfflineChan(geom->Plane(chan.Plane())->View(),
                                  chan.Cell()));
  }

  //......................................................................
  bool AttenCache::Key_t::operator<(const Key_t& rhs) const
  {
    if(fDet < rhs.fDet) return true;
    if(fDet > rhs.fDet) return false;
    return fChan < rhs.fChan;
  }

  //......................................................................
  void AttenCache::Reinit()
  {
    fInitializedData             = false;
    fInitializedMC               = false;
    fInitializedUncalibratedMask = false;
    
    // delete the pointers one by one from the data and MC maps
    for(auto it : fMapData) delete it.second;
    for(auto it : fMapMC  ) delete it.second;

    fMapData.clear();
    fMapMC.clear();
    fMapUncalibratedMask.clear();
    fDiblocks.clear();
  }

} // end namespace calib
////////////////////////////////////////////////////////////////////////