Skip to content
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,15 @@ build/
*.tgz
*.pdf
*.png
*.eig
*~
logs/
jupyter_env
.history
Doc
Inputs/Atmospherics/CAFs
Inputs/DUNE_CAF_files
Inputs/DUNE_spline_files
nohup.out
*.txt
*.log
Expand Down
1 change: 1 addition & 0 deletions Apps/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ foreach(app
SigmaVariation
LikelihoodScan
Variations
ConvertToEigen
)

add_executable(${app} ${app}.cpp)
Expand Down
27 changes: 27 additions & 0 deletions Apps/ConvertToEigen.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
#include "Samples/MaCh3DUNEFactory.h"
#include "Samples/SampleHandlerAtm.h"
#include "Fitters/MaCh3Factory.h"

int main(int argc, char * argv[]) {

MaCh3Utils::MaCh3Usage(argc, argv);
auto fitMan = MaCh3ManagerFactory(argc, argv);

//###############################################################################################################################
//Create SampleHandlerFD objects

ParameterHandlerGeneric* xsec = nullptr;
std::vector<SampleHandlerFD*> DUNEPdfs;
MakeMaCh3DuneInstance(fitMan, DUNEPdfs, xsec);

for (size_t iSample=0;iSample<DUNEPdfs.size();iSample++) {
SampleHandlerAtm* Sample = dynamic_cast<SampleHandlerAtm*>(DUNEPdfs[iSample]);
if (!Sample) {
MACH3LOG_ERROR("Can only convert SampleHandlerAtm CAFs to Eigen matrix binary input file");
throw MaCh3Exception(__FILE__, __LINE__);
}

std::string FileName = "AtmSample_"+Sample->GetTitle()+".eig";
Sample->TransferToEigen(FileName);
}
}
20 changes: 19 additions & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,24 @@ if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD)
set(CMAKE_CXX_STANDARD ${ROOT_CXX_STANDARD})
endif()

#=================== EIGEN

find_package(OpenSSL REQUIRED)

#=================== EIGEN

CPMAddPackage(
NAME Eigen
VERSION 3.4.0
URL https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
DOWNLOAD_ONLY YES
)

if(Eigen_ADDED)
add_library(Eigen INTERFACE IMPORTED)
target_include_directories(Eigen INTERFACE ${Eigen_SOURCE_DIR})
endif()

#=================== DUNEAnaObj
if(NOT DEFINED DUNE_ANAOBJ_BRANCH)
set(DUNE_ANAOBJ_BRANCH v03_06_00)
Expand Down Expand Up @@ -125,7 +143,7 @@ endif()
#If MaCh3 was sourced find it, otherwise use CPM
SET(MaCh3_FOUND FALSE)
find_package(MaCh3 2.2.3 EXACT QUIET)
set(MaCh3_CORE_BRANCH v2.2.3 CACHE STRING "Specify the MaCh3 core branch to use")
set(MaCh3_CORE_BRANCH v2.2.3_patches CACHE STRING "Specify the MaCh3 core branch to use")
message(STATUS "Using MaCh3_CORE_BRANCH: ${MaCh3_CORE_BRANCH}")
option(MaCh3DUNE_Atmos "Whether working with atmospherics or beam oscillations" OFF)

Expand Down
2 changes: 2 additions & 0 deletions Configs/Samples/AtmSample_nueselec.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ Binning:
SampleOptions:
IsELike: yes
ExposureScaling: 1.0
EigenFile: "" #"Inputs/Atmospherics/AtmSample_nueselec.eig"
EigenFileMD5Sum: "" #"964c4c23a6dc52523896c7dc0e264850"

InputFiles:
mtupleprefix: "Inputs/Atmospherics/CAFs/atm_hd"
Expand Down
2 changes: 2 additions & 0 deletions Configs/Samples/AtmSample_numuselec.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ Binning:
SampleOptions:
IsELike: no
ExposureScaling: 1.0
EigenFile: "" #"Inputs/Atmospherics/AtmSample_numuselec.eig"
EigenFileMD5Sum: "" #"da5e519c5565ee315230ce5f71df377d"

InputFiles:
mtupleprefix: "Inputs/Atmospherics/CAFs/atm_hd"
Expand Down
3 changes: 2 additions & 1 deletion Samples/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ set(HEADERS
SampleHandlerBeamND.h
SampleHandlerBeamNDGAr.h
MaCh3DUNEFactory.h
DUNEUtils.h
)
set(SOURCE_FILES
SampleHandlerAtm.cpp
Expand All @@ -25,7 +26,7 @@ if(NOT CPU_ONLY)
set_target_properties(SamplesDUNE PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
endif()

target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE duneanaobj::all MaCh3::All MaCh3DUNECompilerOptions)
target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE duneanaobj::all MaCh3::All Eigen OpenSSL::SSL MaCh3DUNECompilerOptions)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ROOT has TMD5, might be easier than assuming every target has OpenSSL, though they probably do....

target_link_libraries(SamplesDUNE PRIVATE DUNEMaCh3Warnings)

target_include_directories(SamplesDUNE PUBLIC
Expand Down
116 changes: 116 additions & 0 deletions Samples/DUNEUtils.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
#pragma once

#include <iomanip>
#include <vector>
#include <iostream>

#include <Eigen/Dense>
#include <openssl/evp.h>

template<class Matrix>
inline void WriteEigenMatrixToFile(const std::string& filename, const Matrix& matrix){
std::ofstream out(filename, std::ios::out | std::ios::binary | std::ios::trunc);
if(out.is_open()) {
typename Matrix::Index rows=matrix.rows(), cols=matrix.cols();
out.write(reinterpret_cast<char*>(&rows), sizeof(typename Matrix::Index));
out.write(reinterpret_cast<char*>(&cols), sizeof(typename Matrix::Index));
out.write(reinterpret_cast<const char*>(matrix.data()), rows*cols*static_cast<typename Matrix::Index>(sizeof(typename Matrix::Scalar)) );
out.close();
}
else {
MACH3LOG_ERROR("Can not write to file:{}",filename);
throw MaCh3Exception(__FILE__, __LINE__);
}
}

template<class Matrix>
inline void ReadEigenMatrixFromFile(const std::string& filename, Matrix& matrix){
std::ifstream in(filename, std::ios::in | std::ios::binary);
if (in.is_open()) {
typename Matrix::Index rows=0, cols=0;
in.read(reinterpret_cast<char*>(&rows),sizeof(typename Matrix::Index));
in.read(reinterpret_cast<char*>(&cols),sizeof(typename Matrix::Index));
matrix.resize(rows, cols);
in.read(reinterpret_cast<char*>(matrix.data()), rows*cols*static_cast<typename Matrix::Index>(sizeof(typename Matrix::Scalar)) );
in.close();
}
else {
MACH3LOG_ERROR("Can not open binary matrix file:{}",filename);
throw MaCh3Exception(__FILE__, __LINE__);
}
}

// Compute MD5 digest of buffer 'data'. Returns true on success and fills 'digest' (16 bytes).
inline void ComputeMD5(const unsigned char* data, size_t len, unsigned char digest[EVP_MAX_MD_SIZE], unsigned int &digest_len) {
EVP_MD_CTX *ctx = nullptr;

ctx = EVP_MD_CTX_new();
if (!ctx) {
std::cerr << "EVP_MD_CTX_new failed\n";
throw MaCh3Exception(__FILE__, __LINE__);
}

// Use EVP_md5() — still available under OpenSSL 3.0's default provider.
if (1 != EVP_DigestInit_ex(ctx, EVP_md5(), nullptr)) {
std::cerr << "EVP_DigestInit_ex failed\n";
throw MaCh3Exception(__FILE__, __LINE__);
}

if (len > 0) {
if (1 != EVP_DigestUpdate(ctx, data, len)) {
std::cerr << "EVP_DigestUpdate failed\n";
throw MaCh3Exception(__FILE__, __LINE__);
}
}

if (1 != EVP_DigestFinal_ex(ctx, digest, &digest_len)) {
std::cerr << "EVP_DigestFinal_ex failed\n";
throw MaCh3Exception(__FILE__, __LINE__);
}

EVP_MD_CTX_free(ctx);
}

// Read entire file into vector<char>, returns false on error
inline bool ReadFileIntoVector(const std::string &path, std::vector<unsigned char> &out) {
std::ifstream ifs(path, std::ios::binary | std::ios::ate);
if (!ifs) return false;
std::streamsize size = ifs.tellg();
ifs.seekg(0, std::ios::beg);
out.resize(static_cast<size_t>(size));
if (!ifs.read(reinterpret_cast<char*>(out.data()), size)) return false;
return true;
}

// Convert binary digest to lowercase hex string
inline std::string ToHex(const unsigned char* data, size_t len) {
std::ostringstream oss;
oss << std::hex << std::setfill('0');
for (size_t i = 0; i < len; ++i) {
oss << std::setw(2) << static_cast<int>(data[i]);
}
return oss.str();
}

inline bool CompareMD5Sum(std::string KnownSum, std::string FilePathToCheck) {
MACH3LOG_INFO("Determining MD5 Sum for File:{}",FilePathToCheck);

std::vector<unsigned char> Data;
if (!ReadFileIntoVector(FilePathToCheck,Data)) {
throw MaCh3Exception(__FILE__, __LINE__);
}

unsigned char Digest[EVP_MAX_MD_SIZE];
unsigned int Digest_len = 0;
ComputeMD5(Data.data(), Data.size(), Digest, Digest_len);
std::string CalculatedMD5Sum = ToHex(Digest, Digest_len);

MACH3LOG_INFO("Comparing to known MD5 Sum: {}",KnownSum);
if (KnownSum != CalculatedMD5Sum) {
MACH3LOG_ERROR("MD5 Sums are different! Calculated MD5: {}",CalculatedMD5Sum);
throw MaCh3Exception(__FILE__, __LINE__);
}

MACH3LOG_INFO("MD5 Sums are identical!");
return true;
}
78 changes: 64 additions & 14 deletions Samples/SampleHandlerAtm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,9 @@
#include "duneanaobj/StandardRecord/StandardRecord.h"
#pragma GCC diagnostic pop

#include <Eigen/Dense>
#include "DUNEUtils.h"

SampleHandlerAtm::SampleHandlerAtm(std::string mc_version_, ParameterHandlerGeneric* xsec_cov_, const std::shared_ptr<OscillationHandler>& Oscillator_) : SampleHandlerFD(mc_version_, xsec_cov_, Oscillator_) {
KinematicParameters = &KinematicParametersDUNE;
ReversedKinematicParameters = &ReversedKinematicParametersDUNE;
Expand All @@ -17,10 +20,15 @@ SampleHandlerAtm::~SampleHandlerAtm() {
}

void SampleHandlerAtm::Init() {
dunemcSamples.resize(nSamples,dunemc_atm());

IsELike = Get<bool>(SampleManager->raw()["SampleOptions"]["IsELike"],__FILE__,__LINE__);
ExposureScaling = Get<double>(SampleManager->raw()["SampleOptions"]["ExposureScaling"],__FILE__,__LINE__);

if (SampleManager->raw()["SampleOptions"]["EigenFile"]) {
EigenInputFile = Get<std::string>(SampleManager->raw()["SampleOptions"]["EigenFile"],__FILE__,__LINE__);
EigenInputFileMD5Sum = Get<std::string>(SampleManager->raw()["SampleOptions"]["EigenFileMD5Sum"],__FILE__,__LINE__);
} else {
EigenInputFile = "";
}
}

void SampleHandlerAtm::SetupSplines() {
Expand All @@ -33,11 +41,64 @@ void SampleHandlerAtm::SetupWeightPointers() {
MCSamples[i].total_weight_pointers.push_back(MCSamples[i].osc_w_pointer);
MCSamples[i].total_weight_pointers.push_back(&(MCSamples[i].xsec_w));
MCSamples[i].total_weight_pointers.push_back(&(ExposureScaling));
}
}

int SampleHandlerAtm::ReadFromEigen() {
Eigen::MatrixXd Matrix;

MACH3LOG_INFO("Reading Eigen matrix from:{}",EigenInputFile);
CompareMD5Sum(EigenInputFileMD5Sum,EigenInputFile);
ReadEigenMatrixFromFile(EigenInputFile.c_str(),Matrix);

int nEntries = static_cast<int>(Matrix.rows());
dunemcSamples.resize(nEntries);

for (size_t iEvent=0;iEvent<dunemcSamples.size();++iEvent) {
dunemcSamples[iEvent].Target = static_cast<int>(Matrix(iEvent,Variables::Target));
dunemcSamples[iEvent].nupdg = static_cast<int>(Matrix(iEvent,Variables::nupdg));
dunemcSamples[iEvent].nupdgUnosc = static_cast<int>(Matrix(iEvent,Variables::nupdgUnosc));
dunemcSamples[iEvent].rw_isCC = static_cast<int>(Matrix(iEvent,Variables::rw_isCC));
dunemcSamples[iEvent].OscChannelIndex = Matrix(iEvent,Variables::OscChannelIndex);
dunemcSamples[iEvent].rw_erec = Matrix(iEvent,Variables::rw_erec);
dunemcSamples[iEvent].rw_etru = Matrix(iEvent,Variables::rw_etru);
dunemcSamples[iEvent].flux_w = Matrix(iEvent,Variables::flux_w);
dunemcSamples[iEvent].mode = Matrix(iEvent,Variables::mode);
dunemcSamples[iEvent].rw_theta = Matrix(iEvent,Variables::rw_theta);
dunemcSamples[iEvent].rw_truecz = Matrix(iEvent,Variables::rw_truecz);
}


return nEntries;
}

void SampleHandlerAtm::TransferToEigen(std::string FileName) {
Eigen::MatrixXd Matrix = Eigen::MatrixXd(dunemcSamples.size(),Variables::nVariables);

for (size_t iEvent=0;iEvent<dunemcSamples.size();++iEvent) {
Matrix(iEvent,Variables::Target) = dunemcSamples[iEvent].Target;
Matrix(iEvent,Variables::nupdg) = dunemcSamples[iEvent].nupdg;
Matrix(iEvent,Variables::nupdgUnosc) = dunemcSamples[iEvent].nupdgUnosc;
Matrix(iEvent,Variables::rw_isCC) = dunemcSamples[iEvent].rw_isCC;
Matrix(iEvent,Variables::OscChannelIndex) = dunemcSamples[iEvent].OscChannelIndex;
Matrix(iEvent,Variables::rw_erec) = dunemcSamples[iEvent].rw_erec;
Matrix(iEvent,Variables::rw_etru) = dunemcSamples[iEvent].rw_etru;
Matrix(iEvent,Variables::flux_w) = dunemcSamples[iEvent].flux_w;
Matrix(iEvent,Variables::mode) = dunemcSamples[iEvent].mode;
Matrix(iEvent,Variables::rw_theta) = dunemcSamples[iEvent].rw_theta;
Matrix(iEvent,Variables::rw_truecz) = dunemcSamples[iEvent].rw_truecz;
}

MACH3LOG_INFO("Writing Eigen matrix to:{}",FileName);
WriteEigenMatrixToFile(FileName.c_str(),Matrix);
}

int SampleHandlerAtm::SetupExperimentMC() {
// If the Eigen input file is define, read from that. Otherwise read from CAF file
if (EigenInputFile != "") {
int nEntries = ReadFromEigen();
return nEntries;
}

int CurrErrorLevel = gErrorIgnoreLevel;
gErrorIgnoreLevel = kFatal;

Expand Down Expand Up @@ -175,14 +236,3 @@ double SampleHandlerAtm::ReturnKinematicParameter(int KinematicVariable, int iEv
double SampleHandlerAtm::ReturnKinematicParameter(std::string KinematicParameter, int iEvent) {
return *GetPointerToKinematicParameter(KinematicParameter, iEvent);
}

std::vector<double> SampleHandlerAtm::ReturnKinematicParameterBinning(std::string KinematicParameterStr) {
KinematicTypes KinPar = static_cast<KinematicTypes>(ReturnKinematicParameterFromString(KinematicParameterStr));
return ReturnKinematicParameterBinning(KinPar);
}

std::vector<double> SampleHandlerAtm::ReturnKinematicParameterBinning(KinematicTypes KinPar) {
(void)KinPar;
std::vector<double> ReturnVec;
return ReturnVec;
}
Loading