-
Whitney Armstrong authorednew file: tests/elastic_test.cxx modified : tests / elastic_test.sh
933b1ff5
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);
}