PhotonLibrary.h

Go to the documentation of this file.

#ifndef PHOTONLIBRARY_H
#define PHOTONLIBRARY_H

#include "larsim/PhotonPropagation/IPhotonLibrary.h"

#include "TTree.h"
#include "TF1.h"
#include "larsim/Simulation/PhotonVoxels.h"

#include "lardataobj/Utilities/LazyVector.h"


namespace phot{
  
  class PhotonLibrary : public IPhotonLibrary
  {
  public:
    PhotonLibrary() = default;

    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);
    

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

  private:
    
    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;
    
    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)));
    }

    static std::string const OpChannelBranchName;
    
    static size_t ExtractNOpChannels(TTree* tree);

    static int size_t2int(size_t val) {
    return (val <= INT_MAX) ? (int)((ssize_t)val) : -1; }
  };

}

#endif