ShowerCalibrationGaloreFromPID.h

Go to the documentation of this file.

#ifndef LAREXAMPLES_SERVICES_SHOWERCALIBRATIONGALORE_PROVIDERS_SHOWERCALIBRATIONGALOREFROMPID_H
#define LAREXAMPLES_SERVICES_SHOWERCALIBRATIONGALORE_PROVIDERS_SHOWERCALIBRATIONGALOREFROMPID_H


// LArSoft libraries
#include "larexamples/Services/ShowerCalibrationGalore/Providers/ShowerCalibrationGalore.h"

#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"
#include "cetlib_except/exception.h"

// ROOT libraries
#include "Math/Interpolator.h"
#include "TDirectory.h"
#include "TClass.h"
#include "TH1.h"
#include "TObject.h"

// C/C++ standard libraries
#include <string>
#include <vector>
#include <memory> // std::unique_ptr()
#include <type_traits> // std::is_base_of<>
#include <initializer_list>


// forward declarations
class TGraph; // from ROOT


namespace lar {
   namespace example {
      
      // BEGIN ShowerCalibrationGalore -----------------------------------------
      
      namespace details {
         
         template <typename ROOTobj>
         std::unique_ptr<ROOTobj> readROOTobject
            (TDirectory* sourceDir, std::string name);
         
      } // namespace details
      
      
      std::pair<std::string, std::string> splitROOTpath(std::string path);
      
      
      class ShowerCalibrationGaloreFromPID: public ShowerCalibrationGalore {
            public:
         
         //---------------------------------------------------------------------
         struct Config {
            using Name = fhicl::Name;
            using Comment = fhicl::Comment;
            
            fhicl::Atom<std::string> CalibrationFile {
               Name("CalibrationFile"),
               Comment(
                 "path to calibration file and ROOT directory"
                 " (e.g. path/to/file.root:Dir/Dir)"
                 )
            };
            
         }; // struct Config
         
         using parameters_type = fhicl::Table<Config>;
         
         
         //---------------------------------------------------------------------
         ShowerCalibrationGaloreFromPID(Config const& config)
           { readCalibration(config.CalibrationFile()); }
         
         //---------------------------------------------------------------------
         ShowerCalibrationGaloreFromPID(fhicl::ParameterSet const& pset)
           : ShowerCalibrationGaloreFromPID
             (parameters_type(pset, { "service_type", "service_provider" })())
           {}
         
         
         
         //---------------------------------------------------------------------
         virtual float correctionFactor
           (recob::Shower const&, PDGID_t = unknownID) const override;
         
         virtual Correction_t correction
           (recob::Shower const&, PDGID_t = unknownID) const override;
         
         
         virtual std::string report() const override
           { std::ostringstream sstr; reportTo(sstr); return sstr.str(); }
         
         
         void readCalibration(std::string path);
         
         
         template <typename Stream>
         void reportTo(Stream&& out) const;
         
         
         static void verifyOrder(TGraph const* graph);
         
         
         //---------------------------------------------------------------------
            private:
         struct CalibrationInfo_t {
            std::vector<PDGID_t> appliesTo; 
            double minE = -1.; 
            double maxE = -1.; 
            
            std::unique_ptr<ROOT::Math::Interpolator> factor;
            
            std::unique_ptr<ROOT::Math::Interpolator> error;
            
            double evalFactor(double E) const;
            double evalError(double E) const;
            
            bool present() const { return maxE >= 0.; }
            bool uniform() const { return minE == maxE; }
            
            CalibrationInfo_t& applyTo(PDGID_t id);
            CalibrationInfo_t& applyTo(std::initializer_list<PDGID_t> ids);
            
            template <typename Stream>
            void reportTo(Stream&& out) const;
         }; // CalibrationInfo_t
         
         
         CalibrationInfo_t Calibration_pi0; 
         CalibrationInfo_t Calibration_photon; 
         CalibrationInfo_t Calibration_electron; 
         CalibrationInfo_t Calibration_muon; 
         CalibrationInfo_t Calibration_other; 
         
         
         CalibrationInfo_t const& selectCorrection(PDGID_t id) const;
         
         CalibrationInfo_t readParticleCalibration
           (TDirectory* SourceDir, std::string GraphName) const;
            
         CalibrationInfo_t readParticleCalibration
           (TDirectory* SourceDir, std::string GraphName, PDGID_t id) const;
            
         CalibrationInfo_t readParticleCalibration(
            TDirectory* SourceDir, std::string GraphName,
            std::initializer_list<PDGID_t> ids
            ) const;
         
         static TDirectory* OpenROOTdirectory(std::string path);
         
         static std::unique_ptr<ROOT::Math::Interpolator>
         createInterpolator(unsigned int N, double const* x, double const* y);
      
      }; // class ShowerCalibrationGaloreFromPID
      
      // END ShowerCalibrationGalore -----------------------------------------------
      
      
   } // namespace example
} // namespace lar



//------------------------------------------------------------------------------
//--- template implementation
//--- 
template <typename ROOTobj>
std::unique_ptr<ROOTobj> lar::example::details::readROOTobject
   (TDirectory* sourceDir, std::string name)
{
   if (!sourceDir) {
      throw cet::exception("readROOTobject")
        << "Invalid source ROOT directory\n";
   }
   
   // read the object and immediately claim its ownership
   TObject* pObj = sourceDir->Get(name.c_str());
   if (!pObj) {
      throw cet::exception("readROOTobject")
        << "No object '" << name << "' in ROOT directory '"
        << sourceDir->GetPath() << "'\n";
   }

   if (std::is_base_of<TH1, std::decay_t<ROOTobj>>::value)
      static_cast<TH1*>(pObj)->SetDirectory(nullptr);
//   if (std::is_base_of<TTree, std::decay_t<ROOTobj>>::value)
//      static_cast<TTree*>(pObj)->SetDirectory(nullptr);
   std::unique_ptr<TObject> obj(pObj);
   
   // use ROOT to investigate whether this is the right type
   if (!obj->InheritsFrom(ROOTobj::Class())) {
      throw cet::exception("readROOTobject")
        << "Object '" << name << "' in ROOT directory '"
        << sourceDir->GetPath() << "' is a " << obj->IsA()->GetName()
        << ", not derived from " << ROOTobj::Class()->GetName() << "\n";
   }
   
   // now that we are sure obj is of the right type, we transfer ownership
   // from obj to the return value (that here is implicitly constructed)
   return std::unique_ptr<ROOTobj>(static_cast<ROOTobj*>(obj.release()));
} // lar::example::details::readROOTobject()

//------------------------------------------------------------------------------
//--- lar::example::ShowerCalibrationGaloreFromPID
//--- 
template <typename Stream>
void lar::example::ShowerCalibrationGaloreFromPID::reportTo(Stream&& out) const
{
   out << "Corrections for:";
   out << "\n  - neutral pion:      ";
   Calibration_pi0.reportTo(std::forward<Stream>(out));
   out << "\n  - photon:            ";
   Calibration_photon.reportTo(std::forward<Stream>(out));
   out << "\n  - electron/positron: ";
   Calibration_electron.reportTo(std::forward<Stream>(out));
   out << "\n  - muon/antimuon:     ";
   Calibration_muon.reportTo(std::forward<Stream>(out));
   out << "\n  - other (default):   ";
   Calibration_other.reportTo(std::forward<Stream>(out));
   out << "\n";
} // lar::example::ShowerCalibrationGaloreFromPID::report()


//------------------------------------------------------------------------------
template <typename Stream>
void lar::example::ShowerCalibrationGaloreFromPID::CalibrationInfo_t::reportTo
   (Stream&& out) const
{
   if (!present()) {
      out << "not present";
      return;
   }
   if (uniform()) {
      out << "uniform correction " << factor->Eval(minE) << " +/- "
         << error->Eval(minE) << " for all energies";
   }
   else {
      out << "correction valid from E=" << minE
         << " GeV (" << factor->Eval(minE) << " +/- " << error->Eval(minE)
         << ") to E=" << maxE
         << " GeV (" << factor->Eval(maxE) << " +/- " << error->Eval(maxE)
         << ")";
   }
   if (!appliesTo.empty()) {
      out << "; covers particles ID={";
      for (auto id: appliesTo) out << " " << id;
      out << " }";
   }
   
} // lar::example::ShowerCalibrationGaloreFromPID::CalibrationInfo_t::report()


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


#endif // LAREXAMPLES_SERVICES_SHOWERCALIBRATIONGALORE_PROVIDERS_SHOWERCALIBRATIONGALOREFROMPID_H