LinearEnergyAlg.h

Go to the documentation of this file.

#ifndef LARRECO_CALORIMETRY_LINEARENERGYALG_H
#define LARRECO_CALORIMETRY_LINEARENERGYALG_H


// LArSoft libraries
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataalg/DetectorInfo/DetectorProperties.h"
#include "lardataalg/DetectorInfo/DetectorClocks.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" // util::kModBoxA ...

// infrastructure and utilities
// #include "cetlib/exception.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/Assns.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"
#include "fhiclcpp/ParameterSet.h"

// C/C++ standard libraries
#include <vector>
#include <utility> // std::forward()
#include <type_traits> // std::is_same, std::decay_t

namespace calo {
  
  class LinearEnergyAlg {

    public:
    
    struct ModelName {
      static const std::string ModBox;
      static const std::string Birks;
      static const std::string Constant;
    };
    
    enum RecombinationModel_t {
      kModBox,
      kBirks,
      kConstant,
      nRecombinationModel
    };
    
    
    struct RecombinationConfig {
      
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;
      
      bool modelIsBirks() const
        { return Model() == ModelName::Birks; }
      bool modelIsModBox() const
        { return Model() == ModelName::ModBox; }
      bool modelIsConstant() const
        { return Model() == ModelName::Constant; }
      
      
      fhicl::Atom<std::string> Model {
        Name("Model"),
        Comment(std::string("Which recombination model to use: "
          + ModelName::ModBox + ", "
          + ModelName::Birks + ", "
          + ModelName::Constant + ".").c_str()
          )
        };
      
      fhicl::Atom<double> A {
        Name("A"),
        Comment("Parameter \"A\" of box model."),
        fhicl::use_if(this, &RecombinationConfig::modelIsModBox),
        util::kModBoxA
        };

      fhicl::Atom<double> B {
        Name("B"),
        Comment("Parameter \"B\" of box model [kV/cm*(g/cm^2)/MeV]."),
        fhicl::use_if(this, &RecombinationConfig::modelIsModBox),
        util::kModBoxB
        };
      
      fhicl::Atom<double> A3t {
        Name("A3t"),
        Comment("Recombination parameter \"A\" of Birks model."),
        fhicl::use_if(this, &RecombinationConfig::modelIsBirks),
        util::kRecombA
        };

      fhicl::Atom<double> k3t {
        Name("k3t"),
        Comment("Recombination parameter \"k\" of Birks model [kV/cm*(g/cm^2)/MeV]."),
        fhicl::use_if(this, &RecombinationConfig::modelIsBirks),
        util::kRecombk
        };
      
      fhicl::Atom<double> factor {
        Name("factor"),
        Comment("Constant recombination factor for \"constant\" model."),
        fhicl::use_if(this, &RecombinationConfig::modelIsConstant)
        };
      
    }; // RecombinationConfig
    
    
    struct Config {

      using Name = fhicl::Name;
      using Comment = fhicl::Comment;

      fhicl::Atom<bool> UseArea {
        Name("UseArea"),
        Comment("Whether to use the area of hit to count the deposited charges.")
      };

      fhicl::Table<RecombinationConfig> Recombination {
        Name("Recombination"),
        Comment("Parameters of the recombination model")
        };
      
    }; // Config
    
    

    LinearEnergyAlg(Config const& config);
    
    
    LinearEnergyAlg(fhicl::ParameterSet const& pset)
      : LinearEnergyAlg(fhicl::Table<Config>(pset, {})())
      {}



    void setup( detinfo::DetectorProperties const& detproperty, detinfo::DetectorClocks const& detclock, geo::GeometryCore const& geometry )
      { detp = &detproperty; detc = &detclock; geom = &geometry; initialize(); }
    
    template <typename Stream>
    void DumpConfiguration
      (Stream&& out, std::string indent, std::string firstIndent) const;
    
    template <typename Stream>
    void DumpConfiguration(Stream&& out, std::string indent = "") const
      { DumpConfiguration(std::forward<Stream>(out), indent, indent); }

    
    
    template <typename BeginHitIter, typename EndHitIter>
    double CalculateClusterEnergy
      (recob::Cluster const& cluster, BeginHitIter beginHit, EndHitIter endHit)
      const;
    
    template <typename Hits>
    double CalculateClusterEnergy
      (recob::Cluster const& cluster, Hits&& hits) const
      {
        using std::begin; using std::end;
        return CalculateClusterEnergy(cluster, begin(hits), end(hits));
      }
    
    
    std::vector<double> CalculateEnergy(
      std::vector<art::Ptr<recob::Cluster>> const& clusters,
      art::Assns<recob::Cluster, recob::Hit> const& hitsPerCluster
      ) const;
    
    
  private:
    
    struct ModBoxParameters {
      double A = util::kModBoxA;
      double B = util::kModBoxB;
    };
    
    struct BirksParameters {
      double A = util::kRecombA;
      double k = util::kRecombk;
    };
    
    struct ConstantRecombParameters {
      double factor = 1.0; // no sensible default value here
    };
    
    geo::GeometryCore const* geom = nullptr;

    detinfo::DetectorProperties const* detp = nullptr;

    detinfo::DetectorClocks const* detc = nullptr;

    bool   fUseArea;
    int    fRecombModel;
    
    ModBoxParameters recombModBoxParams;
    BirksParameters recombBirksParams;
    ConstantRecombParameters recombConstParams;
    
    // TODO
    // double fRecombA;
    // double fRecombk;
    // double fModBoxA;
    // double fModBoxB;
    // double fRecombFactor;
    double fElectronLifetime;
    double fDeconNorm;

    static constexpr double kWion           = 23.6e-9;  
    // Conversion for energy deposited in GeV to number of ionization electrons produced
    static constexpr double kRecombFactor   = 0.62;     

    void initialize();
    
    double CalculateHitEnergy(recob::Hit const& hit) const;

    double RecombinationCorrection( double dEdx ) const;  

    double ModBoxInverse( double dEdx ) const;

    double BirksInverse( double dEdx ) const;

  };  // class LinearEnergyAlg

}  // namespace calo



//------------------------------------------------------------------------------
//--- template implementation
//---
template <typename BeginHitIter, typename EndHitIter>
double calo::LinearEnergyAlg::CalculateClusterEnergy
  (recob::Cluster const& cluster, BeginHitIter beginHit, EndHitIter endHit)
  const
{
  static_assert( // check the hits type
    std::is_same<std::decay_t<decltype(**beginHit)>, recob::Hit>::value,
    "The hits must be stored as recob::Hit pointers!" // any pointer will do
    );
  
  double E = 0.0; // total cluster energy
  
  for (auto hitIter = beginHit; hitIter != endHit; ++hitIter) {
    
    E += CalculateHitEnergy(**hitIter);
    
  } // for
    
  return E;
  
} // calo::LinearEnergyAlg::CalculateEnergy()


//------------------------------------------------------------------------------
template <typename Stream>
void calo::LinearEnergyAlg::DumpConfiguration
  (Stream&& out, std::string indent, std::string firstIndent) const
{
  
  out << firstIndent << "LinearEnergyAlg configuration:"
    << "\n" << indent << "  use hit " << (fUseArea? "area": "peak amplitude")
      << " for charge estimation"
    << "\n" << indent << "  recombination model: ";
  switch ( fRecombModel ) {
    case kModBox:
      out << ModelName::ModBox
        << "\n" << indent << "    A = " << recombModBoxParams.A
        << "\n" << indent << "    B = " << recombModBoxParams.B;
      break;
    case kBirks:
      out << ModelName::Birks
        << "\n" << indent << "    A = " << recombBirksParams.A
        << "\n" << indent << "    k = " << recombBirksParams.k;
      break;
    case kConstant:
      out << ModelName::Constant
        << "\n" << indent << "    k = " << recombConstParams.factor;
      break;
    default:
      out << "invalid (" << ((int) fRecombModel) << ")!!!";
  } // switch
  out << "\n";
  
} // calo::LinearEnergyAlg::DumpConfiguration()


//------------------------------------------------------------------------------





#endif  // LARRECO_CALORIMETRY_LINEARENERGYALG_H