plot_RFdir_payload.cpp

Go to the documentation of this file.

 
#include <filesystem>
 
#include "util.h"
// CERN ROOT
#include "ROOT/RDataFrame.hxx"
#include "TH2.h"
#include "TVector3.h"
#include "TCanvas.h"
 
// NiceMC
#include "Event.h"
 
int main(int argc, char *argv[]){
 
  std::unique_ptr<TChain> allTree_chain_ptr;
  std::unique_ptr<TChain> passTree_chain_ptr;
  std::filesystem::path PUEO_MC_DATA = argv[1];
  (void) argc; // suppresses unused variable warning
 
  Util::prepare_icefinal_chain(PUEO_MC_DATA, passTree_chain_ptr, allTree_chain_ptr);
 
  TH2D rfdir_hist("weighted passed events",
                  "RF Direction (payload coordinates); #phi; #theta",
                  100, -4,  4,  // horizontal bin number and bounds
                  100, 0, M_PI);
 
  // a lambda expression that fills the above histogram
  auto fill_histogram = [&rfdir_hist](
    const double &path_weight, const double &pos_weight, const double &dir_weight,
    std::vector<nicemc::DetectorEvent> & dE){
 
      double weight = path_weight / (pos_weight * dir_weight);
 
      // detectorEvents[0] is for the main instrument
      TVector3 RFdir_payload = dE[0].RFdir_payload;
 
      // fill the histogram (weighted)
      rfdir_hist.Fill(RFdir_payload.Phi(), RFdir_payload.Theta(), weight);
  };
 
  // use a dataframe to loop through each entry in the passTree
  ROOT::RDataFrame df(*passTree_chain_ptr);
  df.Foreach(
    fill_histogram, 
    {
      "neutrino.path.weight", "loop.positionWeight", "loop.directionWeight",
      "detectorEvents"
    }
  );
 
  TCanvas canvas("", "");
  rfdir_hist.Draw("colz");
  rfdir_hist.SetStats(0);
  canvas.SaveAs("RFdirection_histogram.svg");
 
 
  return 0;
}