effective_area.cpp

#include <regex>
#include "util.h"
#include "effective_area.h"
 
namespace fs = std::filesystem;
 
// constructor, serving as the "main function"
effective_area::effective_area(fs::path p): PUEO_MC_DATA(p) {
  // if there's a run/ with nonexistent IceFinal.root file, have the program exit early.
  if(Util::prepare_icefinal_chain(PUEO_MC_DATA, passTree_chain_ptr, allTree_chain_ptr)) {
    std::cerr << "No further steps will be taken for " << PUEO_MC_DATA << ",\n"
              << "Please fix the missing file problem and retry." << std::endl;
    return;
  } else {
    total_num_neutrino_simulated = allTree_chain_ptr->GetEntries();
    passTree_length = passTree_chain_ptr->GetEntries();
    df = std::make_unique<ROOT::RDataFrame>(*passTree_chain_ptr);
  }
 
  if (passTree_length > 0) {
    effective_area::compute_mean_interaction_length();
 
    if (interaction_length > 0.) {
      effective_area::compute_total_weight();
      effective_area::compute_effective_area();
    }
  }
}
 
void effective_area::compute_mean_interaction_length() {
  auto comp = [this] (const double l){interaction_length += l;};
  df->Foreach(comp, {"interaction.length"});
  interaction_length /= (passTree_length * 1e3 ); // m -> km
}
 
void effective_area::compute_total_weight() {
 
  // a lambda expression to be applied to the DataFrame's weight columns.
  auto count = [this] (const double &w1, const double &w2, const double &w3){
    if (w1 >0 && w2 >0 && w3>0){
      total_num_weighted_passEvents += w1/(w2*w3);
    } else {
      total_num_weighted_passEvents += 0.;
    }
  };
  // loop through all entries in the dataframe to compute the total weight of each entry
  df->Foreach(count,{"neutrino.path.weight", "loop.positionWeight", "loop.directionWeight"});
}
 
void effective_area::compute_effective_area() {
  effective_area_score = ICE_VOLUME * 4 * M_PI * total_num_weighted_passEvents /
                         (total_num_neutrino_simulated * interaction_length);
}