PhotonLibraryHybrid.cxx

Go to the documentation of this file.

#include "larsim/PhotonPropagation/PhotonLibraryHybrid.h"

#include "larsim/Simulation/PhotonVoxels.h"

#include "larcore/Geometry/Geometry.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "TFile.h"
#include "TTree.h"
#include "TVectorD.h"

#include <iostream>

namespace {
  bool fatal(const std::string& msg)
  {
    std::cerr << "FATAL: PhotonLibraryHybrid: " << msg << std::endl;
    abort();
  }
}

namespace phot {
  //--------------------------------------------------------------------
  PhotonLibraryHybrid::PhotonLibraryHybrid(const std::string& fname,
                                           const sim::PhotonVoxelDef& voxdef)
    : fVoxDef(voxdef)
  {
    TFile f(fname.c_str());
    !f.IsZombie() || fatal("Could not open PhotonLibrary " + fname);

    for (int opdetIdx = 0; true; ++opdetIdx) {
      const std::string dirname = TString::Format("opdet_%d", opdetIdx).Data();
      TDirectory* dir = (TDirectory*)f.Get(dirname.c_str());
      if (!dir) break; // Ran out of opdets

      OpDetRecord rec;

      TVectorD* fit = (TVectorD*)dir->Get("fit");
      fit || fatal("Didn't find " + dirname + "/fit in " + fname);
      rec.fit = FitFunc((*fit)[0], (*fit)[1]);

      TTree* tr = (TTree*)dir->Get("tr");
      tr || fatal("Didn't find " + dirname + "/tr in " + fname);
      int vox;
      float vis;
      tr->SetBranchAddress("vox", &vox);
      tr->SetBranchAddress("vis", &vis);

      rec.exceptions.reserve(tr->GetEntries());
      for (int i = 0; i < tr->GetEntries(); ++i) {
        tr->GetEntry(i);
        vox < NVoxels() || fatal("Voxel out of range");
        rec.exceptions.emplace_back(vox, vis);
      }
      rec.exceptions.shrink_to_fit(); // In case we don't trust reserve()

      // In case they weren't sorted in the file
      std::sort(rec.exceptions.begin(), rec.exceptions.end());

      fRecords.push_back(rec);
    } // end for opdetIdx

    !fRecords.empty() || fatal("No opdet_*/ directories in " + fname);

    art::ServiceHandle<geo::Geometry const> geom;
    geom->NOpDets() == fRecords.size() || fatal("Number of opdets mismatch");

    fRecords.shrink_to_fit(); // save memory
  }

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

  //--------------------------------------------------------------------
  int PhotonLibraryHybrid::NVoxels() const
  {
    return fVoxDef.GetNVoxels();
  }

  //--------------------------------------------------------------------
  const float* PhotonLibraryHybrid::GetCounts(size_t vox) const
  {
    // The concept of GetCounts() conflicts fairly badly with how this library
    // works. This is probably the best we can do.

    static float* counts = 0;
    if (!counts) counts = new float[NOpChannels()];
    for (int od = 0; od < NOpChannels(); ++od)
      counts[od] = GetCount(vox, od);
    return counts;
  }

  //--------------------------------------------------------------------
  float PhotonLibraryHybrid::GetCount(size_t vox, size_t opchan) const
  {
    int(vox) < NVoxels() || fatal("GetCount(): Voxel out of range");
    int(opchan) < NOpChannels() || fatal("GetCount(): OpChan out of range");

    const OpDetRecord& rec = fRecords[opchan];

    auto it2 = std::lower_bound(rec.exceptions.begin(), rec.exceptions.end(), vox);
    if (it2->vox == vox) return it2->vis;

    // Otherwise, we use the interpolation, which requires a distance

    static art::ServiceHandle<geo::Geometry const> geom;
    const geo::OpDetGeo& opdet = geom->OpDetGeoFromOpDet(opchan);

    const auto voxvec = fVoxDef.GetPhotonVoxel(vox).GetCenter();

    const double dist = opdet.DistanceToPoint(voxvec);

    return rec.fit.Eval(dist);
  }
}