Skip to content
Snippets Groups Projects
elastic_test.cxx 19.48 KiB
#include "nlohmann/json.hpp"
#include <cmath>
#include <iostream>

#include "ROOT/RDataFrame.hxx"
#include "ROOT/RVec.hxx"

#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TStyle.h"
#include "TSystem.h"
R__LOAD_LIBRARY(libMathMore.so)
R__LOAD_LIBRARY(libGenVector.so)

R__LOAD_LIBRARY(libfmt.so)
#include "fmt/core.h"

#include "THStack.h"

#ifdef __cpp_lib_filesystem
#include <filesystem>
namespace fs = std::filesystem;
#else
#include <experimental/filesystem>
namespace fs = std::experimental::filesystem;
#endif

using RDFNode = decltype(ROOT::RDataFrame{0}.Filter(""));
using Histo1DProxy =
    decltype(ROOT::RDataFrame{0}.Histo1D(ROOT::RDF::TH1DModel{"", "", 128u, 0., 0.}, ""));

struct RDFInfo {
  RDFNode&          df;
  const std::string title;
  RDFInfo(RDFNode& df, std::string_view title) : df{df}, title{title} {}
};

constexpr const double M_P     = .938272;
constexpr const double M_P2    = M_P * M_P;
constexpr const double M_pion  = 0.139;
constexpr const double M_pion2 = M_pion * M_pion;
constexpr const double M_e     = .000511;

// =================================================================================
// Cuts
// =================================================================================
std::string goodTrackSHMS = "P.gtr.dp > -10 && P.gtr.dp < 22";
std::string goodTrackHMS  = "H.gtr.dp > -8 && H.gtr.dp < 8";

std::string piCutSHMS = "P.cal.etottracknorm<1.0";
//std::string piCutSHMS = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && P.cal.etottracknorm<1.0";

std::string eCutHMS = "H.cal.etottracknorm > 0.50 && H.cal.etottracknorm < 2. && "
                      "H.cer.npeSum > 1.";

std::string epiCut = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && "
                     "H.cer.npeSum > 1.0 && H.cal.etottracknorm > 0.6 && "
                     "H.cal.etottracknorm < 2.0 && P.cal.etottracknorm<1.0";

using Pvec3D = ROOT::Math::XYZVector;
using Pvec4D = ROOT::Math::PxPyPzMVector;

// =================================================================================
// reconstruction
// =================================================================================
auto p_pion = [](double px, double py, double pz) {
  return Pvec4D{px * 0.996, py * 0.996, pz * 0.996, M_e};
};
auto p_electron = [](double px, double py, double pz) {
  return Pvec4D{px * 0.994, py * 0.994, pz * 0.994, M_e};
};
auto t = [](const double Egamma, Pvec4D& jpsi) {
  Pvec4D beam{0, 0, Egamma, 0};
  return (beam - jpsi).M2();
};

bool root_file_exists(std::string rootfile) {
  bool found_good_file = false;
  if (!gSystem->AccessPathName(rootfile.c_str())) {
    TFile file(rootfile.c_str());
    if (file.IsZombie()) {
      std::cout << rootfile << " is a zombie.\n";
      std::cout
          << " Did your replay finish?  Check that the it is done before running this script.\n";
      return false;
      // return;
    } else {
      std::cout << " using : " << rootfile << "\n";
      return true;
    }
  }
  return false;
}

void elastic_test(int RunNumber = 6012, int nevents = 50000, int prompt = 0, int update = 0,
                        int default_count_goal = 10000, int redo_timing = 0) {

  // ===============================================================================================
  // Initialization
  // ===============================================================================================
  using nlohmann::json;
  json j;
  {
    std::ifstream json_input_file("db2/run_list_coin.json");
    try {
      json_input_file >> j;
    } catch (json::parse_error) {
      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
      std::quick_exit(-127);
    }
  }

  auto runnum_str = std::to_string(RunNumber);
  if (j.find(runnum_str) == j.end()) {
    std::cout << "Run " << RunNumber << " not found in db2/run_list_coin.json\n";
    std::cout << "Check that run number and replay exists. \n";
    std::cout << "If problem persists please contact Sylvester (217-848-0565)\n";
  }
  double P0_shms_setting = j[runnum_str]["spectrometers"]["shms_momentum"].get<double>();
  double P0_shms         = std::abs(P0_shms_setting);

  bool found_good_file = false;

  std::string rootfile =
      fmt::format("full_online/coin_replay_production_{}_{}.root", RunNumber, nevents);
  found_good_file = root_file_exists(rootfile.c_str());
  if (!found_good_file) {
    rootfile =
        fmt::format("ROOTfiles_volatile/coin_replay_production_{}_{}.root", RunNumber, nevents);
    found_good_file = root_file_exists(rootfile.c_str());
  }
  if (!found_good_file) {
    rootfile = fmt::format("ROOTfiles_csv/coin_replay_production_{}_{}.root", RunNumber, nevents);
    found_good_file = root_file_exists(rootfile.c_str());
  }
  if (!found_good_file) {
    rootfile = fmt::format("ROOTfiles/coin_replay_production_{}_{}.root", RunNumber, nevents);
    found_good_file = root_file_exists(rootfile.c_str());
  }
  if (!found_good_file) {
    std::cout << " Error: suitable root file not found\n";
    return;
  }

  // ===============================================================================================
  // Dataframe
  // ===============================================================================================

  ROOT::EnableImplicitMT(24);

  //---------------------------------------------------------------------------
  // Detector tree
  ROOT::RDataFrame d("T", rootfile);

  //// SHMS Scaler tree
  //ROOT::RDataFrame d_sh("TSP", rootfile);
  //// int N_scaler_events = *(d_sh.Count());

  auto d_coin = d.Filter("fEvtHdr.fEvtType == 4");

  // Good track cuts
  auto dHMSGoodTrack  = d_coin.Filter(goodTrackHMS);
  auto dSHMSGoodTrack = d_coin.Filter(goodTrackSHMS);
  auto dCOINGoodTrack = dHMSGoodTrack.Filter(goodTrackSHMS)
                            .Define("p_electron", p_electron, {"H.gtr.px", "H.gtr.py", "H.gtr.pz"})
                            .Define("p_pion", p_pion, {"P.gtr.px", "P.gtr.py", "P.gtr.pz"});
  // PID cuts
  auto dHMSEl  = dHMSGoodTrack.Filter(eCutHMS);
  auto dSHMSEl = dSHMSGoodTrack.Filter(piCutSHMS);
  auto dCOINEl = dCOINGoodTrack.Filter(eCutHMS + " && " + piCutSHMS);
                     //.Filter(
                     //    [=](double npe, double dp) {
                     //      double p_track = P0_shms * (100.0 + dp) / 100.0;
                     //      // no cerenkov cut needed when momentum is below 2.8 GeV/c
                     //      if (p_track < 2.8) {
                     //        return true;
                     //      }
                     //      return npe > 1.0;
                     //    },
                     //    {"P.hgcer.npeSum", "P.gtr.dp"});

  // Timing cuts
  // Find the timing peak
  // Find the coin peak
  double coin_peak_center = 0;
  if (redo_timing) {
    auto h_coin_time =
        dCOINEl.Histo1D({"coin_time", "coin_time", 8000, 0, 1000}, "CTime.ePositronCoinTime_ROC2");
    h_coin_time->DrawClone();
    int coin_peak_bin = h_coin_time->GetMaximumBin();
    coin_peak_center  = h_coin_time->GetBinCenter(coin_peak_bin);
    std::cout << "COINCIDENCE time peak found at: " << coin_peak_center << std::endl;
  } else {
    //coin_peak_center = 43.0; // run 7240-7241
    coin_peak_center = 23.0; // run 6012
    std::cout << "COINCIDENCE time peak: using pre-calculated value at: " << coin_peak_center
              << std::endl;
    ;
  }
  // timing cut lambdas
  // TODO: evaluate timing cut and offset for random background
  auto timing_cut = [=](double coin_time) { return std::abs(coin_time - coin_peak_center) < 2.; };
  // anti-timing set to 5x width of regular
  auto anti_timing_cut = [=](double coin_time) {
    return std::abs(coin_time - coin_peak_center - 28.) < 10.;
  };
  // timing counts
  auto dHMSElInTime  = dHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
  auto dHMSElRandom  = dHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});
  auto dSHMSElInTime = dSHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
  auto dSHMSElRandom = dSHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});

  auto dCOINElInTime = dCOINEl.Filter(timing_cut, {"CTime.ePiCoinTime_ROC2"});
  auto dCOINElRandom = dCOINEl.Filter(anti_timing_cut, {"CTime.ePiCoinTime_ROC2"});

  // Output root file
  //auto out_file =
  //    new TFile(fmt::format("monitoring/{}/good_csv_counter.root", RunNumber).c_str(), "UPDATE");
  //out_file->cd();

  // =========================================================================================
  // Histograms
  // =========================================================================================
  // 2D correlations
  auto hTheta2DNoCuts = d_coin.Histo2D(
      {"theta2D", "No cuts;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
      "P.gtr.th", "H.gtr.th");
  auto hTheta2DTracking = dCOINGoodTrack.Histo2D(
      {"theta2D", "Cuts: tracking;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
      "P.gtr.th", "H.gtr.th");
  auto hTheta2DPID =
      dCOINEl.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1,
                       .1, 50, -.1, .1},
                      "P.gtr.th", "H.gtr.th");
  auto hTheta2DTiming =
      dCOINElInTime.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50,
                             -.1, .1, 50, -.1, .1},
                            "P.gtr.th", "H.gtr.th");
  // timing
  auto hCoinTimeNoCuts =
      d_coin.Histo1D({"coin_time.NoCuts", "No Cuts;coin_time;counts", 8000, 0, 1000},
                     "CTime.ePositronCoinTime_ROC2");
  auto hCoinTimeTracking = dCOINGoodTrack.Histo1D(
      {"coin_time.Tracking", "Cuts: Tracking;coin_time;counts", 8000, 0, 1000},
      "CTime.ePositronCoinTime_ROC2");
  auto hCoinTimePID =
      dCOINEl.Histo1D({"coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
                      "CTime.ePositronCoinTime_ROC2");
  auto hCoinTimeTiming = dCOINElInTime.Histo1D(
      {"coin_time.Timing", "Cuts: Tracking+PID+Timing;coin_time;counts", 8000, 0, 1000},
      "CTime.ePositronCoinTime_ROC2");

  auto hRandCoinTimePID = dCOINElRandom.Histo1D(
      {"rand_coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
      "CTime.ePositronCoinTime_ROC2");

  // P.gtr.dp
  auto hPdpNoCuts =
      d_coin.Histo1D({"P.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
  auto hPdpTracking = dSHMSGoodTrack.Histo1D(
      {"P.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
  auto hPdpPID = dSHMSEl.Histo1D(
      {"P.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
  auto hPdpTiming = dSHMSElInTime.Histo1D(
      {"P.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
      "P.gtr.dp");
  // P.gtr.th
  auto hPthNoCuts = d_coin.Histo1D(
      {"P.gtr.th.NoCuts", "No Cuts;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
  auto hPthTracking = dSHMSGoodTrack.Histo1D(
      {"P.gtr.th.Tracking", "Cuts: Tracking;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
  auto hPthPID = dSHMSEl.Histo1D(
      {"P.gtr.th.PID", "Cuts: Tracking+PID;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
  auto hPthTiming = dSHMSElInTime.Histo1D(
      {"P.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{SHMS};counts", 200, -0.1, 0.1},
      "P.gtr.th");
  // P.gtr.ph
  auto hPphNoCuts =
      d_coin.Histo1D({"P.gtr.ph.NoCuts", "No Cuts;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
  auto hPphTracking = dSHMSGoodTrack.Histo1D(
      {"P.gtr.ph.Tracking", "Cuts: Tracking;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
  auto hPphPID = dSHMSEl.Histo1D(
      {"P.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
  auto hPphTiming = dSHMSElInTime.Histo1D(
      {"P.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{SHMS};counts", 200, -0.1, 0.1},
      "P.gtr.ph");
  // P.gtr.y
  auto hPyNoCuts =
      d_coin.Histo1D({"P.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "P.gtr.y");
  auto hPyTracking = dSHMSGoodTrack.Histo1D(
      {"P.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "P.gtr.y");
  auto hPyPID =
      dSHMSEl.Histo1D({"P.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "P.gtr.y");
  auto hPyTiming = dSHMSElInTime.Histo1D(
      {"P.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "P.gtr.y");
  // P.cal.etottracknorm
  auto hPcalEPNoCuts =
      d_coin.Histo1D({"P.cal.etottracknorm.NoCuts", "No Cuts;SHMS E/P;counts", 200, -.5, 1.5},
                     "P.cal.etottracknorm");
  auto hPcalEPTracking = dSHMSGoodTrack.Histo1D(
      {"P.cal.etottracknorm.Tracking", "Cuts: Tracking;SHMS E/P;counts", 200, -.5, 1.5},
      "P.cal.etottracknorm");
  auto hPcalEPPID = dSHMSEl.Histo1D(
      {"P.cal.etottracknorm.PID", "Cuts: Tracking+PID;SHMS E/P;counts", 200, -.5, 1.5},
      "P.cal.etottracknorm");
  auto hPcalEPAll = dCOINElInTime.Histo1D(
      {"P.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;SHMS E/P;counts", 200, -.5, 1.5},
      "P.cal.etottracknorm");
  // P.ngcer.npeSum
  auto hPcerNpheNoCuts = d_coin.Histo1D(
      {"P.ngcer.npeSum.NoCuts", "No Cuts;SHMS NGC #phe;counts", 200, -5, 76}, "P.ngcer.npeSum");
  auto hPcerNpheTracking = dSHMSGoodTrack.Histo1D(
      {"P.ngcer.npeSum.Tracking", "Cuts: Tracking;SHMS NGC #phe;counts", 200, -5, 76},
      "P.ngcer.npeSum");
  auto hPcerNphePID = dSHMSEl.Histo1D(
      {"P.ngcer.npeSum.PID", "Cuts: Tracking+PID;SHMS NGC #phe;counts", 200, -5, 76},
      "P.ngcer.npeSum");
  auto hPcerNpheAll = dCOINElInTime.Histo1D(
      {"P.ngcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS NGC #phe;counts", 200, -5, 76},
      "P.ngcer.npeSum");
  // P.hgcer.npeSum
  auto hPhgcerNpheNoCuts = d_coin.Histo1D(
      {"P.hgcer.npeSum.NoCuts", "No Cuts;SHMS HGC #phe;counts", 200, -5, 76}, "P.hgcer.npeSum");
  auto hPhgcerNpheTracking = dSHMSGoodTrack.Histo1D(
      {"P.hgcer.npeSum.Tracking", "Cuts: Tracking;SHMS HGC #phe;counts", 200, -5, 76},
      "P.hgcer.npeSum");
  auto hPhgcerNphePID = dSHMSEl.Histo1D(
      {"P.hgcer.npeSum.PID", "Cuts: Tracking+PID;SHMS HGC #phe;counts", 200, -5, 76},
      "P.hgcer.npeSum");
  auto hPhgcerNpheAll = dCOINElInTime.Histo1D(
      {"P.hgcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS HGC #phe;counts", 200, -5, 76},
      "P.hgcer.npeSum");
  // H.cal.etottracknorm
  auto hHcalEPNoCuts =
      d_coin.Histo1D({"H.cal.etottracknorm.NoCuts", "No Cuts;HMS E/P;counts", 200, -.5, 1.5},
                     "H.cal.etottracknorm");
  auto hHcalEPTracking = dHMSGoodTrack.Histo1D(
      {"H.cal.etottracknorm.Tracking", "Cuts: Tracking;HMS E/P;counts", 200, -.5, 1.5},
      "H.cal.etottracknorm");
  auto hHcalEPPID = dHMSEl.Histo1D(
      {"H.cal.etottracknorm.PID", "Cuts: Tracking+PID;HMS E/P;counts", 200, -.5, 1.5},
      "H.cal.etottracknorm");
  auto hHcalEPAll = dCOINElInTime.Histo1D(
      {"H.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;HMS E/P;counts", 200, -.5, 1.5},
      "H.cal.etottracknorm");
  // H.cer.npeSum
  auto hHcerNpheNoCuts = d_coin.Histo1D(
      {"H.cer.npeSum.NoCuts", "No Cuts;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
  auto hHcerNpheTracking = dHMSGoodTrack.Histo1D(
      {"H.cer.npeSum.Tracking", "Cuts: Tracking;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
  auto hHcerNphePID = dHMSEl.Histo1D(
      {"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
      "H.cer.npeSum");
  auto hHcerNpheAll = dCOINElInTime.Histo1D(
      {"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
      "H.cer.npeSum");
  // H.gtr.dp
  auto hHdpNoCuts =
      d_coin.Histo1D({"H.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
  auto hHdpTracking = dHMSGoodTrack.Histo1D(
      {"H.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
  auto hHdpPID = dHMSEl.Histo1D(
      {"H.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
  auto hHdpTiming = dHMSElInTime.Histo1D(
      {"H.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
      "H.gtr.dp");
  // H.gtr.th
  auto hHthNoCuts = d_coin.Histo1D(
      {"H.gtr.th.NoCuts", "No Cuts;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
  auto hHthTracking = dHMSGoodTrack.Histo1D(
      {"H.gtr.th.Tracking", "Cuts: Tracking;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
  auto hHthPID = dHMSEl.Histo1D(
      {"H.gtr.th.PID", "Cuts: Tracking+PID;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
  auto hHthTiming = dHMSElInTime.Histo1D(
      {"H.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{HMS};counts", 200, -0.1, 0.1},
      "H.gtr.th");
  // H.gtr.ph
  auto hHphNoCuts =
      d_coin.Histo1D({"H.gtr.ph.NoCuts", "No Cuts;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
  auto hHphTracking = dHMSGoodTrack.Histo1D(
      {"H.gtr.ph.Tracking", "Cuts: Tracking;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
  auto hHphPID = dHMSEl.Histo1D(
      {"H.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
  auto hHphTiming = dHMSElInTime.Histo1D(
      {"H.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{HMS};counts", 200, -0.1, 0.1},
      "H.gtr.ph");
  // H.gtr.y
  auto hHyNoCuts =
      d_coin.Histo1D({"H.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "H.gtr.y");
  auto hHyTracking = dHMSGoodTrack.Histo1D(
      {"H.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "H.gtr.y");
  auto hHyPID =
      dHMSEl.Histo1D({"H.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "H.gtr.y");
  auto hHyTiming = dHMSElInTime.Histo1D(
      {"H.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "H.gtr.y");

  // scalers
  //auto total_charge        = d_sh.Max("P.BCM4B.scalerChargeCut");
  //auto shms_el_real_scaler = d_sh.Max("P.pEL_REAL.scaler");
  //auto hms_el_real_scaler  = d_sh.Max("P.hEL_REAL.scaler");
  //auto time_1MHz           = d_sh.Max("P.1MHz.scalerTime");
  //auto time_1MHz_cut       = d_sh.Max("P.1MHz.scalerTimeCut");

  auto yield_all = d.Count();
  // 5 timing cut widths worth of random backgrounds
  auto yield_coin          = d_coin.Count();
  auto yield_HMSGoodTrack  = dHMSGoodTrack.Count();
  auto yield_SHMSGoodTrack = dSHMSGoodTrack.Count();
  auto yield_COINGoodTrack = dCOINGoodTrack.Count();
  auto yield_HMSEl         = dHMSEl.Count();
  auto yield_SHMSEl        = dSHMSEl.Count();
  auto yield_COINEl        = dCOINEl.Count();
  auto yield_HMSElInTime   = dHMSElInTime.Count();
  auto yield_HMSElRandom   = dHMSElRandom.Count();
  auto yield_SHMSElInTime  = dSHMSElInTime.Count();
  auto yield_SHMSElRandom  = dSHMSElRandom.Count();
  auto yield_COINElInTime  = dCOINElInTime.Count();
  auto yield_COINElRandom  = dCOINElRandom.Count();
  auto yield_coin_raw      = dCOINElInTime.Count();
  auto yield_coin_random   = dCOINElRandom.Count();


  // -------------------------------------
  // End lazy eval
  // -------------------------------------
  auto n_coin_good  = *yield_coin_raw - *yield_coin_random / 5.;
  auto n_HMSElGood  = *yield_HMSElInTime - *yield_HMSElRandom / 5;
  auto n_SHMSElGood = *yield_SHMSElInTime - *yield_SHMSElRandom / 5;
  auto n_COINElGood = *yield_COINElInTime - *yield_COINElRandom / 5;

  hPdpNoCuts->DrawCopy();
  //std::cout << "  coin COUNTS : " << *(d_coin.Count()) << "\n";
  //std::cout << " yield_HMSEl : " << *(yield_HMSEl) << "\n";
  std::cout << " yield_COINEl : " << *(yield_COINEl) << "\n";
  //std::cout << "  ALL COUNTS : " << *yield_all << "\n";
  //std::cout << " GOOD COUNTS : " << n_COINElGood << "\n";
  //
  if( 4 != (*(yield_COINEl)) ){
    std::exit(-1);
  }
  std::exit(0);
}