Utilities.cxx

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/Utilities.h"

#include "StandardRecord/Proxy/SRProxy.h"
#include "Utilities/func/Stan.h"

#include "TArrayD.h"
#include "TClass.h"
#include "TDirectory.h"
#include "TH2.h"
#include "TH3.h"
#include "TF1.h"
#include "TMatrixD.h"
#include "TObjString.h"
#include "TString.h"
#include "TVector3.h"
#include "TVectorD.h"

#include <cassert>
#include <cmath>
#include <fstream>
#include <unistd.h>
#include "sys/stat.h"

namespace ana
{
  //----------------------------------------------------------------------
  std::string UniqueName()
  {
    static int N = 0;
    return TString::Format("cafanauniq%d", N++).Data();
  }

  //----------------------------------------------------------------------
  DontAddDirectory::DontAddDirectory()
  {
    fBackup = TH1::AddDirectoryStatus();
    TH1::AddDirectory(false);
  }

  //----------------------------------------------------------------------
  DontAddDirectory::~DontAddDirectory()
  {
    TH1::AddDirectory(fBackup);
  }

  //----------------------------------------------------------------------
  IFDHSilent::IFDHSilent()
  {
    const char* s = getenv("IFDH_SILENT");
    fSet = (s && s == std::string("0"));
    if(!fSet) setenv("IFDH_SILENT", "1", 1);
    if(s) std::cout << "IFDH_SILENT=" << s << std::endl;
  }

  //----------------------------------------------------------------------
  IFDHSilent::~IFDHSilent()
  {
    if(!fSet) unsetenv("IFDH_SILENT");
  }

  //----------------------------------------------------------------------
  FloatingExceptionOnNaN::FloatingExceptionOnNaN(bool enable)
  {
    // Don't want any pending FPEs to trigger when we flip exceptions
    // on. Whoever had them off previously had a reason.
    feclearexcept(FE_INVALID);

    fegetexceptflag(&fBackup, FE_INVALID);

    if(enable)
      feenableexcept(FE_INVALID);
    else
      fedisableexcept(FE_INVALID);
  }

  //----------------------------------------------------------------------
  FloatingExceptionOnNaN::~FloatingExceptionOnNaN()
  {
    fesetexceptflag(&fBackup, FE_INVALID);
  }

  //----------------------------------------------------------------------
  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin, int lastBin)
  {
    if (binSets.size() < 2)
      return std::unique_ptr<TMatrixD>(nullptr);

    if (lastBin < 0)
      lastBin = binSets[0]->GetSize() - 1;  // indexed from 0

    int nBins = lastBin - firstBin + 1;  // firstBin and lastBin are inclusive

    std::vector<double> binMeans(nBins);
    for( const auto & binSet : binSets )
    {
      for ( decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++ )
        binMeans[binIdx] += (*binSet)[binIdx];
    }
    for (decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++)
      binMeans[binIdx] /= binSets.size();


    auto covmx = std::make_unique<TMatrixD>(nBins, nBins);

    for( unsigned int hist_idx = 0; hist_idx < binSets.size(); ++hist_idx )
    {
      // first calculate the weighted sum of squares of the deviations
      for( decltype(nBins) i = 0; i < nBins; i++ )
      {
        double xi = (*(binSets[hist_idx]))[i];
        for( decltype(nBins) k = i; k < nBins; k++ )
        {
          double xk = (*(binSets[hist_idx]))[k];
          (*covmx)[i][k] += (xi - binMeans[i]) * (xk - binMeans[k]);
          if (i != k)
            (*covmx)[k][i] = (*covmx)[i][k];  // covariance matrices are always symmetric
        }
      }
    } // for (hist_idx)

    // now divide by N-1 to get sample covariance
    (*covmx) *= 1./(binSets.size()-1);

    return covmx;
  }

  //----------------------------------------------------------------------
  double LogLikelihood(const TH1* eh, const TH1* oh, bool useOverflow)
  {
    assert(eh->GetNbinsX() == oh->GetNbinsX());

    double chi = 0;

    int bufferBins = useOverflow? 2 : 1;

    for(int i = 0; i < eh->GetNbinsX()+bufferBins; ++i){
      const double e = eh->GetBinContent(i);
      const double o = oh->GetBinContent(i);

      chi += LogLikelihood(e, o);
    }

    return chi;
  }

  //----------------------------------------------------------------------
  TH2F* ExpandedHistogram(const std::string& title,
                          int nbinsx, double xmin, double xmax, bool xlog,
                          int nbinsy, double ymin, double ymax, bool ylog)
  {
    DontAddDirectory guard;

    if(xlog){xmin = log(xmin); xmax = log(xmax);}
    if(ylog){ymin = log(ymin); ymax = log(ymax);}

    // How wide the bins will be once we're done
    const double xwidth = (xmax-xmin)/(nbinsx-1);
    const double ywidth = (ymax-ymin)/(nbinsy-1);

    // Move the bin edges so that the limits occur at the centres
    xmin -= xwidth/2; ymin -= ywidth/2;
    xmax += xwidth/2; ymax += ywidth/2;

    std::vector<double> xedges(nbinsx+1);
    std::vector<double> yedges(nbinsy+1);

    for(int i = 0; i <= nbinsx; ++i){
      xedges[i] = xmin + (xmax-xmin)*i/double(nbinsx);
      if(xlog) xedges[i] = exp(xedges[i]);
    }
    for(int i = 0; i <= nbinsy; ++i){
      yedges[i] = ymin + (ymax-ymin)*i/double(nbinsy);
      if(ylog) yedges[i] = exp(yedges[i]);
    }

    return new TH2F(UniqueName().c_str(), title.c_str(),
                    nbinsx, &xedges.front(),
                    nbinsy, &yedges.front());
  }

  //----------------------------------------------------------------------
  std::string FindPackageDir( std::string Dir )
  {
#ifdef NOVACMAKE
    const char* pub = getenv("NOVASOFT_DIR");
    const char* priv = getenv("NOVASOFT_DIR");
#else
    const char* pub = getenv("FW_BASE");
    const char* priv = getenv("FW_RELEASE_BASE");
#endif

    if(priv && priv != std::string(".")){
      const std::string ret = std::string(priv)+"/"+Dir;
      struct stat junk;
      if(stat(ret.c_str(), &junk) == 0) return ret;
    }
    assert(pub);
    return std::string(pub)+"/"+Dir;
  }

  //----------------------------------------------------------------------
  std::string FindCAFAnaDir()
  {
        return FindPackageDir("CAFAna");
  }

  //----------------------------------------------------------------------
  std::vector<std::string> LoadFileList(const std::string& listfile)
  {
    std::vector<std::string> ret;
    
    std::ifstream is(listfile);
    if(!is.good()){
      std::cerr << "Can't open file list '" << listfile << "'. Aborting." << std::endl;
      abort();
    }

    while(!is.eof()){
      std::string fname;
      is >> fname;
      if(!fname.empty()) ret.push_back(fname);
    }
    return ret;
  }

  //----------------------------------------------------------------------
  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir)
  {
    std::map<std::string, std::string> ret;

    TIter next(dir->GetListOfKeys());
    while(TObject* key = next()){
      TObject* obj = dir->Get(key->GetName());
      assert(obj);

      const TString className = obj->ClassName();
      if(className == "TObjString"){
        ret[key->GetName()] = ((TObjString*)obj)->GetString();
      }
      else{
        std::cerr << "Unexpected object " << key->GetName() << " of type " << className << " while looking for metadata. Ignoring" << std::endl;
      }
    }
    
    return ret;
  }

  //----------------------------------------------------------------------
  void CombineMetadata(std::map<std::string, std::string>& base,
                       const std::map<std::string, std::string>& add,
                       std::set<std::string>& mask)
  {
    for(auto it: add){
      const std::string& key = it.first;

      // Needs special handling anyway, leave it blank.
      if(key == "parents") continue;

      // Accumulate the runs list
      if(key == "runs"){
        const std::string& r1 = base[key];
        const std::string& r2 = it.second;

        assert(!r2.empty());

        // when concatenating MC files,
        // it's possible to get duplicates
        // (run+subrun the same, only differing
        //  in cycle number, which isn't part
        //  of this metadata field)
        // don't repeat, because SAM gets angry.
        if (r1.find(r2.substr(1, r2.size() - 2)) != std::string::npos)
          continue;

        if(r1.empty()){
          base[key] = r2;
          continue;
        }

        // "[foo]" + "[bar]"
        std::string sum = r1+&r2[1]; // "[foo]bar]"
        sum[r1.size()-1] = ',';      // "[foo,bar]"
        base[key] = sum;
        continue;
      }

      if(base.find(key) == base.end()){
        // If it's new, add it
        base[key] = it.second;
      }
      else{
        if(key == "simulated.number_of_spills" ||
           key == "event_count" ||
           key == "online.totalevents"){
          // These two fields should be accumulated
          base[key] = TString::Format("%d",
                                      atoi(base[key].c_str()) +
                                      atoi(it.second.c_str())).Data();
        }
        else{
          // If it's a clash, record it
          if(base[key] != it.second) mask.insert(key);
        }
      }
    }
  }


  //----------------------------------------------------------------------
  void WriteCAFMetadata(TDirectory* dir,
                        const std::map<std::string, std::string>& meta)
  {
    TDirectory* tmp = gDirectory;
    dir->cd();

    for(auto it: meta){
            TObjString str(it.second.c_str());
            str.Write(it.first.c_str());
    }

    dir->Save();

    tmp->cd();
  }

  //----------------------------------------------------------------------
  TVector3 NuMIBeamDirection(caf::Det_t det)
  {
    // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
    // CAF files and fitting the resulting distribution to a gaussian.
    // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
    // and containment cuts.
    // CAF Files used:
    //   NearDet:
    //    prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
    //   FarDet:
    //    prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08 
    //   Unknown for NDOS

    switch(det){
    case caf::kFARDET:  return TVector3(-6.833e-05,
                                        +6.388e-02,
                                        +9.980e-1).Unit();
    case caf::kNDOS:    return TVector3(-3.845e-4,
                                        +6.517e-2,
                                        +9.967e-1).Unit();
    case caf::kNEARDET: return TVector3(-8.424e-04,
                                        -6.174e-02,
                                        +9.981e-1).Unit();
    default:
      std::cerr<< "Aborting. Bad geometry configuration: " << det << "\n";
      abort();
    }
  }

  //----------------------------------------------------------------------
  bool RunningOnGrid()
  {
    static bool cache;
    static bool cache_set = false;
    if(!cache_set){
      cache = (getenv("_CONDOR_SCRATCH_DIR") != 0);
      cache_set = true;
    }

    return cache;
  }

  //----------------------------------------------------------------------
  size_t Stride(bool allow_default)
  {
    static int cache = -1;

    if(cache < 0){
      char* env = getenv("CAFANA_STRIDE");
      if(env){
        cache = std::atoi(env);
      }
      else{
        if(allow_default){
          cache = 1;
        }
        else{
          std::cout << "Stride() called, but CAFANA_STRIDE is not set (--stride not passed?)" << std::endl;
          abort();
        }
      }
    }

    return cache;
  }

  //----------------------------------------------------------------------
  size_t Offset(bool allow_default)
  {
    static int cache = -1;

    if(cache < 0){
      char* env = getenv("CAFANA_OFFSET");
      if(env){
        cache = std::atoi(env);
      }
      else{
        if(allow_default){
          cache = 0;
        }
        else{
          std::cout << "Offset() called, but CAFANA_OFFSET is not set (--offset not passed?)" << std::endl;
          abort();
        }
      }
    }

    return cache;
  }

  //----------------------------------------------------------------------
  int Limit()
  {
    static int cache = 0;

    if(cache == 0){
      char* env = getenv("CAFANA_LIMIT");
      if(env){
        cache = std::atoi(env);
      }
      else{
        cache = -1;
      }
    }

    return cache;
  }

  //----------------------------------------------------------------------
  size_t JobNumber()
  {
    if(!RunningOnGrid()){
      std::cout << "JobNumber() called, but we are not running on the grid" << std::endl;
      abort();
    }

    return Offset(false);
  }

  //----------------------------------------------------------------------
  size_t NumJobs()
  {
    if(!RunningOnGrid()){
      std::cout << "NumJobs() called, but we are not running on the grid" << std::endl;
      abort();
    }

    return Stride(false);
  }

  //----------------------------------------------------------------------
  std::size_t Hash(const std::vector<std::size_t> &vals)
  {
    std::size_t seed = 0;

    for (auto v : vals)
      seed ^= std::hash<std::size_t>{}(v) + 0x9e3779b9 + (seed<<6) + (seed>>2);
    return seed;
  }

  //----------------------------------------------------------------------
  std::size_t Hash(const caf::SRSpillProxy *sr)
  {
    return Hash({sr->run, sr->subrun, sr->evt});
  }

  //----------------------------------------------------------------------
  std::size_t Hash(const caf::SRProxy *sr)
  {
    return Hash({Hash(&(sr->spill)), std::size_t(sr->hdr.cycle), sr->hdr.subevt});
  }

  //----------------------------------------------------------------------
  FitToFourier::FitToFourier(TH1* h, double xlo, double xhi, int NOsc)
    : fHist(h), fxlo(xlo), fxhi(xhi), fNOsc(NOsc)
  {
  }

  //----------------------------------------------------------------------
  FitToFourier::~FitToFourier()
  {
  }

  //----------------------------------------------------------------------
  double FitToFourier::operator()(double *x, double *par) const
  {
    double x0 = x[0];
    double val = par[0];
    for (int i = 1; i <= fNOsc; i++)
      val += par[2*i-1]*sin(i*M_PI*x0) + par[2*i]*cos(i*M_PI*x0);
    return val;
  }

  //----------------------------------------------------------------------
  TF1* FitToFourier::Fit() const
  {
    std::vector<double> s(fNOsc, 0.);
    std::vector<double> c(fNOsc, 0.);
    int nBins = 0;
    for(int i = 1; i <= fHist->GetNbinsX(); ++i){
      const double x = M_PI * fHist->GetXaxis()->GetBinCenter(i);
      const double y = fHist->GetBinContent(i);

      if(y == 0) continue;
      ++nBins;

      for(int n = 0; n <= fNOsc; ++n){
        s[n] += y * sin(n*x);
        c[n] += y * cos(n*x);
      }
    }

    for(int n = 0; n <= fNOsc; ++n){
      s[n] *= 2./nBins;
      c[n] *= 2./nBins;
    }

    TF1* f = new TF1(UniqueName().c_str(), this, fxlo, fxhi, 2*fNOsc+1);

    f->SetParameter(0, c[0]/2);
    for(int n = 1; n <= fNOsc; ++n){
      f->SetParameter(n*2-1, s[n]);
      f->SetParameter(n*2,   c[n]);
    }

    // Because ROOT is having problems drawing f if I don't
    double min = fHist->GetMinimum();
    double max = fHist->GetMaximum();
    f->GetYaxis()->SetRangeUser(0.8*min, 1.2*max);
    return f;
  }
}