NOvARwgtInterface.cxx

Go to the documentation of this file.

/*
 * NOvARwgtInterface.cxx:
 *  Create NOvARwgt objects from StandardRecords
 *
 *  Created on: Jan. 24, 2019
 *      Author: J. Wolcott <jwolcott@fnal.gov>
 */

#include "StandardRecord/NOvARwgt/NOvARwgtInterface.h"

#include <unordered_map>
#include <unordered_set>

#include "nugen/NuReweight/ReweightLabels.h"

#include "StandardRecord/StandardRecord.h"
#include "StandardRecord/SRNeutrino.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include "NOvARwgt/rwgt/EventRecord.h"

// use novarwgt namespace because this is effectively a NOvARwgt file
// (just ported here so that whenever we update SR this doesn't get out of date)
namespace novarwgt
{
        // declared here, not in the .h, so that it doesn't escape this file
        namespace
        {
                const std::unordered_map<rwgt::EReweightLabel, novarwgt::ReweightKnob> RW_KNOB_NAME_MAP
                {
                        {rwgt::fReweightMaCCQE,  novarwgt::kKnob_MaCCQE},
                        {rwgt::fReweightMaCCRES, novarwgt::kKnob_MaCCRES},
                };

                // there's some slight skew between the ordering of these in the GENIE & SR versions
                const std::unordered_map<caf::mode_type_, novarwgt::ReactionType> REACTION_ENUM_MAP
                {
                        {caf::kUnknownMode, novarwgt::kScNull},
                        {caf::kQE, novarwgt::kScQuasiElastic},
                        {caf::kRes, novarwgt::kScResonant},
                        {caf::kDIS, novarwgt::kScDeepInelastic},
                        {caf::kCoh, novarwgt::kScCoherent},
                        {caf::kCohElastic, novarwgt::kScCoherentElas},
                        {caf::kElectronScattering, novarwgt::kScNuElectronElastic},
                        {caf::kIMDAnnihilation, novarwgt::kScIMDAnnihilation},

                        {caf::kInverseBetaDecay, novarwgt::kScInverseBetaDecay},
                        {caf::kGlashowResonance, novarwgt::kScGlashowResonance},
                        {caf::kAMNuGamma, novarwgt::kScAMNuGamma},
                        {caf::kMEC, novarwgt::kScMEC},

                        {caf::kDiffractive, novarwgt::kScDiffractive},

                        // these two don't have GENIE versions...
                        //{caf::kEM, novarwgt::kScNull},
                        //{caf::kWeakMix, novarwgt::kScNull },

                };

                const std::unordered_map<caf::generator_, novarwgt::Generator> GENERATOR_ID_MAP
    {
      {caf::kGENIE, novarwgt::kGENIE}
    };

                // End-of-program warn-er
                class UnmatchedKnobWarn
    {
      public:
        UnmatchedKnobWarn() = default;

        ~UnmatchedKnobWarn()
        {
          if (fUnmatchedKnobs.empty())
            return;

          std::cerr << "WARNING: no corresponding NOvARwgt knob found for the following rwgt::ReweightLabel indices (found in events):" << std::endl;
          for (const auto & knob : fUnmatchedKnobs)
            std::cerr << "  " << knob << ",";
          std::cerr << std::endl;
          std::cerr << "Stored weights for them were ignored." << std::endl;
        }

        void AddKnob(std::size_t newKnob)
        {
          fUnmatchedKnobs.insert(newKnob);
        }

      private:
        std::unordered_set<std::size_t> fUnmatchedKnobs;
    };
        }

        /// Copy information out of an SRNeutrino
        template <typename T>
        const novarwgt::EventRecord& ConvertSRTruth(const T * nu, bool forceNoCache)
        {
                static_assert(std::is_same<T,caf::SRNeutrino>::value || std::is_same<T,caf::SRNeutrinoProxy>::value,
                                "novarwgt::ConvertSRTruth() can only be used with SRNeutrino or SRNeutrinoProxy");

    // hold a cache of the last StandardRecord we converted to a novarwgt::EventRecord
    // so that syst operations that use the same event repeatedly
    //  don't convert it lots of times in a row
    static novarwgt::EventRecord sCachedEvt;

    bool matches =  std::abs(nu->E - sCachedEvt.Enu) < 0.0001 && std::abs(nu->q2 - sCachedEvt.Q2()) < 0.0001
                    && std::abs(nu->y - sCachedEvt.y) < 0.0001 && (nu->iscc == sCachedEvt.isCC);

    if (matches && !forceNoCache)
      return sCachedEvt;

    sCachedEvt.Reset();

    sCachedEvt.generator = GENERATOR_ID_MAP.at(nu->generator);
    for (const auto & v : nu->genVersion)
      sCachedEvt.generatorVersion.push_back(v);
                sCachedEvt.generatorConfigStr = nu->genConfigString;

                sCachedEvt.nupdg      = nu->pdg;
                sCachedEvt.isCC       = nu->iscc;
                sCachedEvt.reaction   = REACTION_ENUM_MAP.at(caf::mode_type_(int(nu->mode)));  // nu->mode was originally a caf::mode_type_
                sCachedEvt.struckNucl = nu->hitnuc;

                sCachedEvt.Enu  = nu->E;
                sCachedEvt.q    = TLorentzVector(1, 0, 0, nu->y * nu->E);  // need to make up the direction since we didn't store it
                sCachedEvt.q.SetRho(sqrt(nu->q2 + sCachedEvt.q.E()*sCachedEvt.q.E()));
        sCachedEvt.y    = nu->y;
                sCachedEvt.W    = sqrt(nu->W2);

                sCachedEvt.A = nu->tgtA;

                sCachedEvt.npiplus = nu->npiplus;
                sCachedEvt.npizero = nu->npizero;
                sCachedEvt.npiminus = nu->npiminus;

                for (const auto & part : nu->prim)
    {
                  sCachedEvt.fsPartMult[part.pdg]++;  // default constructor initializes to 0
                  sCachedEvt.fsPartKE[part.pdg] += (1 - 1./part.p.Gamma()) * part.p.E;
    }

                // the following doesn't work because SRMCReweightProxy owns a VectorProxy<SRGenieWeightsProxy>,
                // and a VectorProxy fundamentally contains a vector<T*> instead of a vector<T>.
//    static_assert(sizeof(novarwgt::ReweightVals) == sizeof(caf::SRGenieWeights),
//                  "The NOvARwgt-StandardRecord interface only works if novarwgt::ReweightVals and caf::SRGenieWeights have the same size!");
//              std::vector<novarwgt::ReweightVals> systWgts(nu->rwgt.genie.size());
//              std::memcpy(systWgts.data(), nu->rwgt.genie.data(), nu->rwgt.genie.size() * sizeof(novarwgt::ReweightVals));
//      sCachedEvt.genieWeights = ReweightList(std::move(systWgts));

    static UnmatchedKnobWarn warn;
                // since we have to dereference the elements of the vector to copy them anyway,
                // we may as well just copy them one by one.
                sCachedEvt.genieWeights.resize(novarwgt::kLastKnob);
                for (std::size_t idx = 0; idx < nu->rwgt.genie.size(); idx++)
                {
                  // slightly less typing
                        const auto & stdRecWgts = nu->rwgt.genie[idx];
                        auto rwgtIdx = static_cast<rwgt::ReweightLabel_t>(idx);
                        try
      {
                          auto novaRwIdx =  novarwgt::kNugenKnobTranslationTable.at(rwgtIdx);
        sCachedEvt.genieWeights[novaRwIdx]
            = {stdRecWgts.minus2sigma, stdRecWgts.minus1sigma, stdRecWgts.plus1sigma, stdRecWgts.plus2sigma};
      }
                        catch(std::out_of_range & e)
      {
                          warn.AddKnob(rwgtIdx);
      }
                }

                if (!nu->isvtxcont)
                        sCachedEvt.expectNoWeights = true;

                // GENIE reweight table maker doesn't make weights for nu-on-e events with Q^2 < 0
                // since the GENIE reweight calcs tend to flip out
                if (nu->rwgt.genie.empty() && nu->mode == caf::kElectronScattering)
                  sCachedEvt.expectNoWeights = true;

                // in prod5.1 we introduced strange/charm CCQE production, and the z-expansion knobs barf on those events.
                // so we simultaneously disabled those knobs for those events.
                // unfortunately that meant they're just missing from the end of the reweight vector...
                // so NOvARwgt just sees them as missing, and then it tries to invoke GENIE to calculate them,
                // causing a crash all over again (but at CAFAna time instead of production time).
                // the underlying issue in MakeGENIEReweightTable was subsequently fixed so that it won't happen in future production,
                // but for now we need to work around it and explicitly set the relevant weights to 1.
                // someday when we're no longer using prod5.1 files we should be able to remove this workaround...
                // note the strange/charm baryon test is inside the innermost test just because it's the most expensive.
                //
                // for an unknown reason any QE events on hydrogen experienced the same fate,
                // even though the branches in the MakeGENIEReweightTable code don't appear to induce it.
                // (it's appropriate, however, because QE-on-H has no need for ANY form factors since the nucleon is free.)
                // investigation ongoing; but the workaround is the same for now.
                bool fixZexpKnobs = false;
                if (nu->iscc && nu->mode == caf::kQE && !sCachedEvt.genieWeights.IsSet(novarwgt::kKnob_ZNormCCQE))
    {
                  if (nu->tgtA == 1)
                    fixZexpKnobs = true;
                  else
      {
        for (const auto &part : nu->prefsi)
        {
          if (std::abs(part.pdg) > 3000 && std::abs(part.pdg) < 5000)
          {
            fixZexpKnobs = true;
            break;
          }
        }
      } // else --> nu->tgtA != 1
    } // if (iscc || is QE || !IsSet(kKnob_ZNormCCQE)

    if (fixZexpKnobs)
    {
      const std::vector <novarwgt::ReweightKnob> knobs_to_fix = {
          novarwgt::kKnob_ZNormCCQE,
          novarwgt::kKnob_ZExpA1CCQE,
          novarwgt::kKnob_ZExpA2CCQE,
          novarwgt::kKnob_ZExpA3CCQE,
          novarwgt::kKnob_ZExpA4CCQE,
          novarwgt::kKnob_AxFFCCQEshape
      };

      for (const auto &knob : knobs_to_fix)
        sCachedEvt.genieWeights[knob] = {1.0, 1.0, 1.0, 1.0};
    }

    return sCachedEvt;

        }
        template const novarwgt::EventRecord& ConvertSRTruth(const caf::SRNeutrino *, bool);
        template const novarwgt::EventRecord& ConvertSRTruth(const caf::SRNeutrinoProxy *, bool);

        /// Copy information out of a StandardRecord
        template <typename T>
        const novarwgt::EventRecord& ConvertSREvent(const T * sr, bool forceNoCache)
        {
                static_assert(std::is_same<T,caf::StandardRecord>::value || std::is_same<T,caf::SRProxy>::value,
                              "novarwgt::ConvertSRTruth() can only be used with SRNeutrino or SRNeutrinoProxy");

                static novarwgt::EventRecord emptyRecord;
                if (sr->mc.nnu != 1)
                        return emptyRecord;

                return ConvertSRTruth( &(sr->mc.nu[0]), forceNoCache );


        }
        template const novarwgt::EventRecord& ConvertSREvent(const caf::StandardRecord *, bool);
        template const novarwgt::EventRecord& ConvertSREvent(const caf::SRProxy *, bool);
}