print_genie_weights.C File Reference

#include <iostream>
#include <string>
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"

Go to the source code of this file.

Functions
const std::string	genie_producer_label ("generator")
const std::string	eventweight_producer_label ("genieeventweightTest")
void	print_weights (const std::string &filename)

Function Documentation

const std::string eventweight_producer_label

(

"genieeventweightTest"

)

Referenced by print_weights().

const std::string genie_producer_label

(

"generator"

)

Referenced by print_weights().

void print_weights

(

const std::string &

filename

)

Definition at line 15 of file print_genie_weights.C.

References simb::MCParticle::E(), eventweight_producer_label(), evwgh::MCEventWeight::fWeight, genie_producer_label(), simb::MCTruth::GetNeutrino(), simb::MCNeutrino::Nu(), and w.

{

  std::vector<std::string> filenames{filename};

  gallery::Event ev(filenames);

  size_t event_count = 1;
  for (ev.toBegin(); !ev.atEnd(); ++ev) {

    std::cout << "art event " << event_count << '\n';

    // Use the GENIE producer label here
    auto mctruth_handle = ev.getValidHandle<std::vector<simb::MCTruth>>(genie_producer_label);
    auto gtruth_handle = ev.getValidHandle<std::vector<simb::GTruth>>(genie_producer_label);

    // Use the EventWeight producer label here
    auto weights_handle =
      ev.getValidHandle<std::vector<evwgh::MCEventWeight>>(eventweight_producer_label);

    // Loop through these objects for each neutrino vertex in the event
    for (size_t v = 0u; v < mctruth_handle->size(); ++v) {

      const simb::MCTruth& mc = mctruth_handle->at(v);
      const simb::GTruth& gt = gtruth_handle->at(v);
      const evwgh::MCEventWeight& mc_weights = weights_handle->at(v);

      // Extract desired information
      const auto& nu = mc.GetNeutrino();
      if (nu.Nu().NumberTrajectoryPoints() > 0) {
        double E_nu = nu.Nu().E(0);
        bool cc = (nu.CCNC() == 0);
        if (cc)
          std::cout << "CC";
        else
          std::cout << "NC";

        int mode = nu.Mode();
        if (mode == 0) {
          if (cc)
            std::cout << "QE";
          else
            std::cout << "EL";
        }
        else if (mode == 1)
          std::cout << "RES";
        else if (mode == 2)
          std::cout << "DIS";
        else if (mode == 3)
          std::cout << "COH";
        else if (mode == 10)
          std::cout << "MEC";
        else
          std::cout << " other";
        std::cout << " interaction with neutrino energy = " << E_nu << " GeV\n";
      }

      // Loop over all of the weights in the MCEventWeight object
      std::cout << "Weights\n";
      for (const auto& pair : mc_weights.fWeight) {
        std::string knob_name = pair.first;
        std::vector<double> weights = pair.second;

        std::cout << "  " << knob_name << ":\n";
        for (size_t u = 0u; u < weights.size(); ++u) {
          double w = weights.at(u);
          std::cout << "    universe #" << u << " has weight = " << w << '\n';
        }
      }
    }

    ++event_count;
  }
}