PhotonLibrary.h

Go to the documentation of this file.

#ifndef PHOTONLIBRARY_H
#define PHOTONLIBRARY_H

#include "larsim/PhotonPropagation/IPhotonLibrary.h"

#include "larsim/Simulation/PhotonVoxels.h"

#include "TF1.h"
class TTree;

#include "lardataobj/Utilities/LazyVector.h"

#include <limits> // std::numeric_limits
#include <optional>

namespace art {
  class TFileDirectory;
}

namespace phot {

  class PhotonLibrary : public IPhotonLibrary {
  public:
    PhotonLibrary(art::TFileDirectory* pDir = nullptr);

    TTree* ProduceTTree() const;

    virtual float GetCount(size_t Voxel, size_t OpChannel) const override;
    void SetCount(size_t Voxel, size_t OpChannel, float Count);

    float GetTimingPar(size_t Voxel, size_t OpChannel, size_t parnum) const;
    void SetTimingPar(size_t Voxel, size_t OpChannel, float Count, size_t parnum);

    //    TF1& GetTimingTF1(size_t Voxel, size_t OpChannel) const;
    void SetTimingTF1(size_t Voxel, size_t OpChannel, TF1 func);

    virtual float GetReflCount(size_t Voxel, size_t OpChannel) const override;
    void SetReflCount(size_t Voxel, size_t OpChannel, float Count);

    virtual float GetReflT0(size_t Voxel, size_t OpChannel) const override;
    void SetReflT0(size_t Voxel, size_t OpChannel, float reflT0);

    virtual float const* GetCounts(size_t Voxel) const override;
    const std::vector<float>* GetTimingPars(size_t Voxel) const;
    TF1* GetTimingTF1s(size_t Voxel) const;

    virtual float const* GetReflCounts(size_t Voxel) const override;
    virtual float const* GetReflT0s(size_t Voxel) const override;

    bool hasTiming() const { return fHasTiming; }

    virtual bool hasReflected() const override { return fHasReflected; }

    virtual bool hasReflectedT0() const override { return fHasReflectedT0; }

    void StoreLibraryToFile(std::string LibraryFile,
                            bool storeReflected = false,
                            bool storeReflT0 = false,
                            size_t storeTiming = 0) const;
    void LoadLibraryFromFile(std::string LibraryFile,
                             size_t NVoxels,
                             bool storeReflected = false,
                             bool storeReflT0 = false,
                             size_t storeTiming = 0,
                             int maxrange = 200);
    void CreateEmptyLibrary(size_t NVoxels,
                            size_t NChannels,
                            bool storeReflected = false,
                            bool storeReflT0 = false,
                            size_t storeTiming = 0);

    // --- BEGIN --- Metadata: voxel information -------------------------------

    bool hasVoxelDef() const { return fVoxelDef.has_value(); }

    sim::PhotonVoxelDef const& GetVoxelDef() const
    {
      assert(fVoxelDef);
      return *fVoxelDef;
    }

    void SetVoxelDef(sim::PhotonVoxelDef const& voxelDef) { fVoxelDef = voxelDef; }

    // --- END --- Metadata: voxel information ---------------------------------

  private:
    virtual int NOpChannels() const override { return fNOpChannels; }
    virtual int NVoxels() const override { return fNVoxels; }

    virtual bool isVoxelValid(size_t Voxel) const override { return isVoxelValidImpl(Voxel); }

    bool fHasReflected = false; 
    bool fHasReflectedT0 =
      false; 

    size_t fHasTiming =
      0; 

    // fLookupTable[unchecked_index(Voxel, OpChannel)] = Count
    // for each voxel, all NChannels() channels are stored in sequence
    util::LazyVector<float> fLookupTable;
    util::LazyVector<float> fReflLookupTable;
    util::LazyVector<float> fReflTLookupTable;
    util::LazyVector<std::vector<float>> fTimingParLookupTable;
    util::LazyVector<TF1> fTimingParTF1LookupTable;
    std::string fTimingParFormula;
    size_t fTimingParNParameters;

    size_t fNOpChannels;
    size_t fNVoxels;

    std::optional<sim::PhotonVoxelDef> fVoxelDef;

    art::TFileDirectory* fDir = nullptr;

    bool isVoxelValidImpl(size_t Voxel) const { return Voxel < fNVoxels; }

    size_t uncheckedIndex(size_t Voxel, size_t OpChannel) const
    {
      return Voxel * fNOpChannels + OpChannel;
    }

    float uncheckedAccess(size_t Voxel, size_t OpChannel) const
    {
      return fLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float& uncheckedAccess(size_t Voxel, size_t OpChannel)
    {
      return fLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float uncheckedAccessRefl(size_t Voxel, size_t OpChannel) const
    {
      return fReflLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float& uncheckedAccessRefl(size_t Voxel, size_t OpChannel)
    {
      return fReflLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float uncheckedAccessReflT(size_t Voxel, size_t OpChannel) const
    {
      return fReflTLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float& uncheckedAccessReflT(size_t Voxel, size_t OpChannel)
    {
      return fReflTLookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    float uncheckedAccessTimingPar(size_t Voxel, size_t OpChannel, size_t parnum) const
    {
      return fTimingParLookupTable[uncheckedIndex(Voxel, OpChannel)][parnum];
    }

    float& uncheckedAccessTimingPar(size_t Voxel, size_t OpChannel, size_t parnum)
    {
      return fTimingParLookupTable[uncheckedIndex(Voxel, OpChannel)][parnum];
    }

    TF1& uncheckedAccessTimingTF1(size_t Voxel, size_t OpChannel)
    {
      return fTimingParTF1LookupTable[uncheckedIndex(Voxel, OpChannel)];
    }

    const TF1& uncheckedAccessTimingTF1(size_t Voxel, size_t OpChannel) const
    {
      // note that this will produce a segmentation fault if the formula is not there
      return *(fTimingParTF1LookupTable.data_address(uncheckedIndex(Voxel, OpChannel)));
    }

    void LoadMetadata(TDirectory& srcDir);

    void StoreMetadata() const;

    static std::string const OpChannelBranchName;

    static size_t ExtractNOpChannels(TTree* tree);

    static int size_t2int(size_t val)
    {
      constexpr size_t int_max_as_size_t{std::numeric_limits<int>::max()};
      return val <= int_max_as_size_t ? static_cast<int>(val) : -1;
    }
  };

}

#endif