VarVals.h

Go to the documentation of this file.

/*
 * \file VarVals.h
 *   \brief Base container for the Vars that constitute an event.
 *
 *  Created on: \date June 13, 2016
 *  Original author:  \author J. Wolcott <jwolcott@fnal.gov>
 *
 */

#ifndef COVARIANCEMATRIXFIT_CORE_VARVALS_H_
#define COVARIANCEMATRIXFIT_CORE_VARVALS_H_

#include <string>
#include <cstdint>

#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "CovarianceMatrixFit/dataProducts/Constants.h"
#include "CovarianceMatrixFit/dataProducts/Structs.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"

namespace cmf
{
  struct TruthVars{

    TruthVars()
    : fTrueE              (std::numeric_limits<float>::lowest())
    , fTruePDG            (std::numeric_limits<float>::lowest())
    , fTrueCCNC           (std::numeric_limits<float>::lowest())
    , fTrueIntType        (std::numeric_limits<float>::lowest())
    , fTruePDGOrig        (std::numeric_limits<float>::lowest())
    , fTrueParentPDG      (std::numeric_limits<float>::lowest())
    , fTrueParentPT       (std::numeric_limits<float>::lowest())
    , fTrueParentPZ       (std::numeric_limits<float>::lowest())
    , fTrueParentDecay    (std::numeric_limits<float>::lowest())
    , fTrueParentTargetPDG(std::numeric_limits<float>::lowest())
    , fTrueHitNuc         (std::numeric_limits<float>::lowest())
    , fTrueIntMode        (std::numeric_limits<float>::lowest())
    {}

    float fTrueE;               ///< True nu energy
    float fTruePDG;             ///< true PDG
    float fTrueCCNC;            ///< true cc vs nc
    float fTrueIntType;         ///< true interaction type
    float fTruePDGOrig;         ///< true PDG of original neutrino
    float fTrueParentPDG;       ///< true PDG of neutrino parent
    float fTrueParentPT;        ///< true p_T of neutrino parent
    float fTrueParentPZ;        ///< true p_Z of neutrino parent
    float fTrueParentDecay;     ///< true parent decay mode
    float fTrueParentTargetPDG; ///< true parent pdg code off the target
    float fTrueHitNuc;          ///< true hit nucleus
    float fTrueIntMode;         ///< true interaction mode
  };

  class EventId{

  public:
    EventId()
    : run   (std::numeric_limits<int>::lowest())
    , subrun(std::numeric_limits<int>::lowest())
    , event (std::numeric_limits<int>::lowest())
    , cycle (std::numeric_limits<int>::lowest())
    , slice (std::numeric_limits<int>::lowest())
    {}

    EventId(int r,
            int s,
            int e,
            int c,
            int slc)
    : run   (r)
    , subrun(s)
    , event (e)
    , cycle (c)
    , slice (slc)
    {}

    int run;
    int subrun;
    int event;
    int cycle;
    int slice;

    friend std::ostream& operator << (std::ostream & o, cmf::EventId const& evid);

    bool operator<(EventId const& evid)  const;
    bool operator==(EventId const& evid) const;

  };

  struct DataVars{

    // initialize the Fake_Weight to be 1, it will be reset if Fake_Data are generated
    DataVars()
    : fNuE_CVN         (std::numeric_limits<float>::lowest())
    , fNuE_NumMichel   (std::numeric_limits<float>::lowest())
    , fNu_RecoE        (std::numeric_limits<float>::lowest()) 
    , fHad_RecoE       (std::numeric_limits<float>::lowest())
    , fLep_RecoE       (std::numeric_limits<float>::lowest())
    , fLep_RecoE_MCFrac(0.)
    , fRecoQ2          (std::numeric_limits<float>::lowest())
    , fFake_Weight     (1.)
    {}

    DataVars(float nuecvn,
             float nuenummichel,
             float nurece,
             float hadrece,
             float leprece,
             float leprecemcfrac,
             float recoQ2,
             float fakewgt)
    : fNuE_CVN         (nuecvn)
    , fNuE_NumMichel   (nuenummichel)
    , fNu_RecoE        (nurece)
    , fHad_RecoE       (hadrece)
    , fLep_RecoE       (leprece)
    , fLep_RecoE_MCFrac(leprecemcfrac)
    , fRecoQ2          (recoQ2)
    , fFake_Weight     (fakewgt)
    {}

    float fNuE_CVN;          ///<
    float fNuE_NumMichel;    ///<
    float fNu_RecoE;         ///<
    float fHad_RecoE;        ///<
    float fLep_RecoE;        ///<
    float fLep_RecoE_MCFrac; ///< fraction of leptonic energy in muon catcher
    float fRecoQ2;           ///< reconstructed Q^2
    float fFake_Weight;      ///< Weight for fake data events
  };

  // set the default values to be unity, ie no weighting applied
  struct WeightVars{

    WeightVars()
    : fXSecCVPPFX_Weight            (1.)
    , fRPACCQE_Weight               (1.)
    , fRPARES_Weight                (1.)
    , fDISvnCC1pi_Weight            (1.)
    , fEmpiricalMEC_Weight          (1.)
    , fEmpiricalMECtoGENIEQE_Weight (1.)
    , fEmpiricalMECtoGENIERES_Weight(1.)
    , fOtherDIS_Weight              (1.)
    , fXSecCV2020_Weight            (1.)
    , fPPFXFluxCV_Weight            (1.)
    {}

    float fXSecCVPPFX_Weight;             ///< Deprecated, Was Tufts weight for SA
    float fRPACCQE_Weight;                ///< Deprecated
    float fRPARES_Weight;                 ///< To be used for systematic evaluation ONLY
    float fDISvnCC1pi_Weight;             ///< Deprecated
    float fEmpiricalMEC_Weight;           ///< Deprecated
    float fEmpiricalMECtoGENIEQE_Weight;  ///< Deprecated
    float fEmpiricalMECtoGENIERES_Weight; ///< Deprecated
    float fOtherDIS_Weight;               ///< Deprecated
    float fXSecCV2020_Weight;             ///< For 2020 Ana, store weights seperately
    float fPPFXFluxCV_Weight;             ///< For 2020 Ana, store weights seperately
  };

  // struct to hold weights of systematic parameters with shifts related to
  // different sigma shifts in the underlying parameter
  struct SystVar {

    SystVar(uint8_t type)
    : fType(type)
    , fMinus2Sigma(1.)
    , fMinus1Sigma(1.)
    , fPlus1Sigma(1.)
    , fPlus2Sigma(1.)
    {}

    SystVar(uint8_t type,
            float minus2,
            float minus1,
            float plus1,
            float plus2)
    : fType(type)
    , fMinus2Sigma(minus2)
    , fMinus1Sigma(minus1)
    , fPlus1Sigma(plus1)
    , fPlus2Sigma(plus2)
    {}

    uint8_t             Type() const { return fType; }
    std::vector<float>  SigmasAsVec() const;
    float               SigmaWeight(uint8_t const& sigmaLevel) const;

    bool operator==(cmf::SystVar const& other) const { return fType == other.Type(); }
    bool operator <(cmf::SystVar const& other) const { return fType  < other.Type(); }
    
    friend std::ostream& operator << (std::ostream & o, cmf::SystVar const& sv);

    uint8_t fType;        ///< key for this parameter
    float   fMinus2Sigma; ///< weight associated with this change in the underlying parameter
    float   fMinus1Sigma; ///< weight associated with this change in the underlying parameter
    float   fPlus1Sigma;  ///< weight associated with this change in the underlying parameter
    float   fPlus2Sigma;  ///< weight associated with this change in the underlying parameter
  };

  //----------------------------------------------------------------------------
  inline std::vector<float> SystVar::SigmasAsVec() const
  {
    return std::move(std::vector<float>({fMinus2Sigma,
                                         fMinus1Sigma,
                                         fPlus1Sigma,
                                         fPlus2Sigma}));
  }

  //----------------------------------------------------------------------------
  inline float SystVar::SigmaWeight(uint8_t const& sigmaLevel) const
  {
         if(sigmaLevel == cmf::kMinus2Sigma) return fMinus2Sigma;
    else if(sigmaLevel == cmf::kMinus1Sigma) return fMinus1Sigma;
    else if(sigmaLevel == cmf::kPlus1Sigma ) return fPlus1Sigma;
    else if(sigmaLevel == cmf::kPlus2Sigma ) return fPlus2Sigma;

    // default to no change in weight
    return 1.;
  }

  typedef std::vector<cmf::SystVar> SystVarColl;
  

  /// \brief Base container for the MC related Vars that constitute an event.
  class MCVarVals
  {
  private:
    cmf::TruthVars   fTruthVars;
    cmf::WeightVars  fWeightVars;
    cmf::SystVarColl fSystematicShifts;

  public:
    MCVarVals()
    {}

    MCVarVals(cmf::TruthVars  const& tv,
              cmf::WeightVars const& wv)
    : fTruthVars (tv)
    , fWeightVars(wv)
    {}
    
    MCVarVals(cmf::TruthVars   const& tv,
              cmf::WeightVars  const& wv,
              cmf::SystVarColl const& gv)
    : fTruthVars (tv)
    , fWeightVars(wv)
    , fSystematicShifts(gv)
    {}


    /// Access a Var by name.
    float val_at(uint8_t const& varkey) const;

    std::map<float, float> SystVarWgts(uint8_t const& varkey) const;

    friend std::ostream& operator << (std::ostream & o, cmf::MCVarVals const& mcvv);

    // The following functions should only be used for making eventlist trees
    cmf::TruthVars   const& TruthVars()         const { return fTruthVars;        }
    cmf::WeightVars  const& WeightVars()        const { return fWeightVars;       }
    cmf::SystVarColl const& SystematicShifts()  const { return fSystematicShifts; }

  };

  /// \brief Base container for the MC related Vars that constitute an event.
  class DataVarVals
  {
  private:

    float NuRecoEVal(cmf::MetaData const& md) const;

    cmf::DataVars fDataVars;

  public:
    DataVarVals () = default;

    DataVarVals(cmf::DataVars const& dv)
    : fDataVars(dv)
    {}

    /// Access a Var by name.
    /// Throws if the var doesn't exist
    float val_at(uint8_t       const& varkey,
                 cmf::MetaData const& md) const;

    void  set_val_at(uint8_t const& varkey,
                     float   const& val);

    friend std::ostream& operator << (std::ostream & o, cmf::DataVarVals const& dvv);

    // The following function should only be used for making eventlist trees
    cmf::DataVars const& DataVars() const { return fDataVars; }

  };

  static float RecoEnergy(float                leptonicEnergy,
                          float                hadronicEnergy,
                          cmf::MetaData const& md);

  static float NuMuEnergy(float leptonicEnergy,
                          float hadronicEnergy);

  static float NuEEnergy(float                  leptonicEnergy,
                         float                  hadronicEnergy,
                         cmf::BeamType_t const& beamType);

  static float NCEnergy(float                  leptonicEnergy,
                        float                  hadronicEnergy,
                        cmf::BeamType_t const& beamType,
                        cmf::DetType_t  const& detector);

} /* namespace cmf */

//----------------------------------------------------------------------------
inline float cmf::RecoEnergy(float                leptonicEnergy,
                             float                hadronicEnergy,
                             cmf::MetaData const& md)
{
  if     (md.IsNuMuSelected()) return cmf::NuMuEnergy(leptonicEnergy, hadronicEnergy               );
  else if(md.IsNCSelected()  ) return cmf::NCEnergy  (leptonicEnergy, hadronicEnergy, md.BeamType(), md.detector);
  else if(md.IsNuESelected() ) return cmf::NuEEnergy (leptonicEnergy, hadronicEnergy, md.BeamType());
  else
    throw cet::exception("DataVarVals")
        << "Unknown selection type "
        << md.selectionType
        << " requested to be used in calculating reconstructed energy";

  return std::numeric_limits<float>::lowest();
}

//----------------------------------------------------------------------------
inline float cmf::NuMuEnergy(float leptonicEnergy,
                             float hadronicEnergy)
{
  return leptonicEnergy + hadronicEnergy;

}

//----------------------------------------------------------------------------
inline float cmf::NuEEnergy(float                  leptonicEnergy,
                            float                  hadronicEnergy,
                            cmf::BeamType_t const& beamType)
{
  // p0 and p3 are 0.0 for both FHC and RHC
  // 2020
  float p0 = 0.0;
  float p1 = (beamType == kFHC) ? 1.01777 : 0.988258;
  float p2 = (beamType == kFHC) ? 1.10868 : 1.20084;
  float p3 = 0.0;
  float p4 = (beamType == kFHC) ? 1.43541e-03 : 1.92904e-07;
  float p5 = (beamType == kFHC) ? 1.09628e-01 : 1.20704e-07;
  float fr = (beamType == kFHC) ? 1. / 1.0355622 : 1. / 1.0111873;

  return fr *(( p5 * hadronicEnergy * hadronicEnergy) +
              ( p4 * leptonicEnergy * leptonicEnergy) +
              ( p3 * leptonicEnergy * hadronicEnergy) +
              ( p2 * hadronicEnergy)                  +
              ( p1 * leptonicEnergy)                  +
              p0 );
}

//----------------------------------------------------------------------------
inline float cmf::NCEnergy(float                      leptonicEnergy,
                           float                      hadronicEnergy,
                           cmf::BeamType_t     const& beamType      ,
                           cmf::DetType_t      const& detector      )
{
  // Note: leptonicEnergy must be the EM-like energy, as calculated by kEME_2020
  // and hadronicEnergy must be the hadronic-like energy, as calculated by kHADE_2020
  // If we don't have a tuned correction return the total calorimetric energy
  if (beamType != kFHC ||
      (detector != cmf::kNEARDET && detector != cmf::kFARDET))
    return leptonicEnergy + hadronicEnergy;

  // Set up correction parameters
  double pars[4][6] = { {1.049, 0.795, 0.8409, 0.17,  0.82, -1.00},  // ND FHC
                        {1.025, 0.797, 0.9162, 0.53, -0.26, -1.48},  // FD FHC
                        {1.000, 1.000, 1.0000, 0.00,  0.00,  0.00},  // ND RHC
                        {1.000, 1.000, 1.0000, 0.00,  0.00,  0.00} };// FD RHC

  // set detector/beam index
  int detCur = (detector == cmf::kFARDET) + ((beamType == kRHC)<<1); // bitshift magic

  // return the energy estimator
  return (leptonicEnergy / pars[detCur][0] + hadronicEnergy / pars[detCur][1])
    / (pars[detCur][2] + pars[detCur][3] * std::pow(leptonicEnergy + hadronicEnergy + pars[detCur][4],2) * std::exp(pars[detCur][5] * (leptonicEnergy + hadronicEnergy)));
}
#endif /* COVARIANCEMATRIXFIT_CORE_VARVALS_H_ */