feat: move main to MaCh3 core v2.5.0

.gitignore

@@ -11,4 +11,6 @@ jupyter_env
 Doc
 Inputs/Atmospherics/CAFs
 nohup.out
-output*.txt
+*.txt
+*.log
+.cache/

Apps/EventRates.cpp

@@ -14,62 +14,60 @@
 
 #include "Samples/MaCh3DUNEFactory.h"
 #include "Samples/StructsDUNE.h"
+#include "Fitters/MaCh3Factory.h"
 
-void Write1DHistogramsToFile(std::string OutFileName, std::vector<TH1*> Histograms) {
+void Write1DHistogramsToFile(std::string OutFileName, std::vector<std::unique_ptr<TH1>>& Histograms) {
   auto OutputFile = std::unique_ptr<TFile>(TFile::Open(OutFileName.c_str(), "RECREATE"));
   OutputFile->cd();
-  for(auto Hist : Histograms){
+  for(auto& Hist : Histograms){
     Hist->Write();
   }
   OutputFile->Close();
 }
 
-void Write1DHistogramsToPdf(std::string OutFileName, std::vector<TH1*> Histograms) {
+void Write1DHistogramsToPdf(std::string OutFileName, std::vector<std::unique_ptr<TH1>>& Histograms) {
   //Remove root from end of file
   OutFileName.erase(OutFileName.find('.'));
   OutFileName+=".pdf";
 
   auto c1 = std::unique_ptr<TCanvas>(new TCanvas("c1", "c1", 800, 600));
   c1->cd();
   c1->Print(std::string(OutFileName+"[").c_str());
-  for(auto Hist : Histograms){
+  for(auto& Hist : Histograms){
     Hist->Draw("HIST");
     c1->Print(OutFileName.c_str());
   }
   c1->Print(std::string(OutFileName+"]").c_str());
 }
 
 int main(int argc, char * argv[]) {
-  MaCh3Utils::MaCh3Usage(argc, argv);
-  auto fitMan = std::unique_ptr<manager>(new manager(argv[1]));
+  M3::Utils::MaCh3Usage(argc, argv);
+  auto FitManager = MaCh3ManagerFactory(argc, argv);
 
   //###############################################################################################################################
-  //Create SampleHandlerFD objects
-
-  ParameterHandlerGeneric* xsec = nullptr;
-  ParameterHandlerOsc* osc = nullptr;
-
-  std::vector<SampleHandlerFD*> DUNEPdfs;
-  MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc);
+  //Create SampleHandlerBase objects
+  auto [param_handler, samples] = MaCh3DuneFactory(FitManager);
 
   //###############################################################################################################################
   //Perform reweight and print total integral
 
-  std::vector<TH1*> DUNEHists;
-  for(auto Sample : DUNEPdfs){
-    Sample->Reweight();
-    if (Sample->GetNDim() == 1)
-      DUNEHists.push_back(Sample->Get1DHist());
-    else if (Sample->GetNDim() == 2)
-      DUNEHists.push_back(Sample->Get2DHist());
-
-    std::string EventRateString = fmt::format("{:.2f}", Sample->Get1DHist()->Integral());
-    MACH3LOG_INFO("Event rate for {} : {:<5}", Sample->GetTitle(), EventRateString);
+  std::vector<std::unique_ptr<TH1>> DUNEHists;
+  for(auto& handler : samples){
+    if (!handler){
+      MACH3LOG_ERROR("Sample not set up correctly");
+      throw MaCh3Exception(__FILE__, __LINE__);
+    }
+    handler->Reweight();
+    for (int iSample=0; iSample<handler->GetNSamples(); iSample++) {
+      DUNEHists.push_back(M3::Clone(handler->GetMCHist(iSample)));
 
-    Sample->PrintIntegral();
+      std::string EventRateString = fmt::format("{:.2f}", handler->GetMCHist(iSample)->Integral());
+      MACH3LOG_INFO("Event rate for {} : {:<5}", handler->GetSampleTitle(iSample), EventRateString);
+      handler->PrintIntegral(iSample);
+    }
   }
 
-  std::string OutFileName = GetFromManager<std::string>(fitMan->raw()["General"]["OutputFile"], "EventRatesOutput.root");
+  std::string OutFileName = GetFromManager<std::string>(FitManager->raw()["General"]["OutputFile"], "EventRatesOutput.root");
   Write1DHistogramsToFile(OutFileName, DUNEHists); 
   Write1DHistogramsToPdf(OutFileName, DUNEHists);
 
@@ -80,24 +78,26 @@ int main(int argc, char * argv[]) {
   MACH3LOG_INFO("========================================================================");
   MACH3LOG_INFO("Oscillation Mode Breakdown:");
 
-  for(auto Sample : DUNEPdfs) {
-    MACH3LOG_INFO("======================");
-    int nOscChannels = Sample->GetNMCSamples();
-    for (int iOscChan=0;iOscChan<nOscChannels;iOscChan++) {
-      std::vector< KinematicCut > SelectionVec;
-
-      KinematicCut SelecChannel;
-      SelecChannel.ParamToCutOnIt = Sample->ReturnKinematicParameterFromString("OscillationChannel");
-      SelecChannel.LowerBound = iOscChan;
-      SelecChannel.UpperBound = iOscChan+1;
-      SelectionVec.push_back(SelecChannel);
-
-      TH1* Hist = Sample->Get1DVarHist(Sample->GetXBinVarName(),SelectionVec);
-      MACH3LOG_INFO("{:<20} : {:<20} : {:<20.2f}",Sample->GetTitle(),Sample->GetFlavourName(iOscChan),Hist->Integral());
+  for(auto handler : samples) {
+    for (int iSample = 0; iSample < handler->GetNSamples(); iSample++) {
+      MACH3LOG_INFO("======================");
+      int nOscChannels = handler->GetNOscChannels(iSample);
+      for (int iOscChan=0;iOscChan<nOscChannels;iOscChan++) {
+        std::vector< KinematicCut > SelectionVec;
+
+        KinematicCut SelecChannel;
+        SelecChannel.ParamToCutOnIt = handler->ReturnKinematicParameterFromString("OscillationChannel");
+        SelecChannel.LowerBound = iOscChan;
+        SelecChannel.UpperBound = iOscChan+1;
+        SelectionVec.push_back(SelecChannel);
+
+        auto Hist = handler->Get1DVarHist(iSample, handler->GetKinVarName(iSample, 0),SelectionVec);
+        MACH3LOG_INFO("{:<20} : {:<20} : {:<20.2f}",handler->GetSampleTitle(iSample),handler->GetFlavourName(iSample, iOscChan),Hist->Integral());
+      }
+
+      auto Hist = handler->Get1DVarHist(iSample, handler->GetKinVarName(iSample, 0));
+      MACH3LOG_INFO("{:<20} : {:<20.2f}",handler->GetSampleTitle(iSample),Hist->Integral());
     }
-
-    TH1* Hist = Sample->Get1DVarHist(Sample->GetXBinVarName());
-    MACH3LOG_INFO("{:<20} : {:<20.2f}",Sample->GetTitle(),Hist->Integral());
   }
 
   //###############################################################################################################################
@@ -107,26 +107,28 @@ int main(int argc, char * argv[]) {
   MACH3LOG_INFO("========================================================================");
   MACH3LOG_INFO("Interaction Mode Breakdown:");
 
-  for(auto Sample : DUNEPdfs) {
-    MACH3LOG_INFO("======================");
+  for(auto handler : samples) {
+    for (int iSample = 0; iSample < handler->GetNSamples(); iSample++) {
+      MACH3LOG_INFO("======================");
 
-    MaCh3Modes* Modes = Sample->GetMaCh3Modes();
-    int nModeChannels = Modes->GetNModes();
-    for (int iModeChan=0;iModeChan<nModeChannels;iModeChan++) {
-      std::vector< KinematicCut > SelectionVec;
+      MaCh3Modes* Modes = handler->GetMaCh3Modes();
+      int nModeChannels = Modes->GetNModes();
+      for (int iModeChan=0;iModeChan<nModeChannels;iModeChan++) {
+        std::vector< KinematicCut > SelectionVec;
 
-      KinematicCut SelecChannel;
-      SelecChannel.ParamToCutOnIt = Sample->ReturnKinematicParameterFromString("Mode");
-      SelecChannel.LowerBound = iModeChan;
-      SelecChannel.UpperBound = iModeChan+1;
-      SelectionVec.push_back(SelecChannel);
+        KinematicCut SelecChannel;
+        SelecChannel.ParamToCutOnIt = handler->ReturnKinematicParameterFromString("Mode");
+        SelecChannel.LowerBound = iModeChan;
+        SelecChannel.UpperBound = iModeChan+1;
+        SelectionVec.push_back(SelecChannel);
 
-      TH1* Hist = Sample->Get1DVarHist(Sample->GetXBinVarName(),SelectionVec);
-      MACH3LOG_INFO("{:<20} : {:<20} : {:<20.2f}",Sample->GetTitle(),Modes->GetMaCh3ModeName(iModeChan),Hist->Integral());
-    }
+        auto Hist = handler->Get1DVarHist(iSample, handler->GetKinVarName(iSample,0), SelectionVec);
+        MACH3LOG_INFO("{:<20} : {:<20} : {:<20.2f}",handler->GetSampleTitle(iSample),Modes->GetMaCh3ModeName(iModeChan),Hist->Integral());
+      }
 
-    TH1* Hist = Sample->Get1DVarHist(Sample->GetXBinVarName());
-    MACH3LOG_INFO("{:<20} : {:<20.2f}",Sample->GetTitle(),Hist->Integral());
+      auto Hist = handler->Get1DVarHist(iSample, handler->GetKinVarName(iSample,0));
+      MACH3LOG_INFO("{:<20} : {:<20.2f}",handler->GetSampleTitle(iSample),Hist->Integral());
+    }
   }
 
   //###############################################################################################################################

Apps/Fit.cpp

@@ -12,28 +12,17 @@
 #include <TColor.h>
 #include <TMath.h>
 
-#include "Fitters//mcmc.h"
+#include "Fitters/MaCh3Factory.h"
 #include "Samples/MaCh3DUNEFactory.h"
 
 int main(int argc, char * argv[]) {
 
-  // ----------------------- OPTIONS ---------------------------------------- //
-  if(argc == 1){
-    MACH3LOG_INFO("Usage: bin/jointFitDUNEBeam Configs/config.yaml");
-    return 1;
-  }
-
-  manager *FitManager = new manager(argv[1]);
+  auto FitManager = MaCh3ManagerFactory(argc, argv);
   auto OutputFileName = FitManager->raw()["General"]["OutputFile"].as<std::string>();
 
-  ParameterHandlerGeneric* xsec = nullptr;
-  ParameterHandlerOsc* osc = nullptr;
-
   //####################################################################################
-  //Create samplePDFSKBase Objs
-
-  std::vector<SampleHandlerFD*> DUNEPdfs;
-  MakeMaCh3DuneInstance(FitManager, DUNEPdfs, xsec, osc); 
+  //Create sample handler + parameter_handler objects
+  auto [param_handler, samples] = MaCh3DuneFactory(FitManager);
 
   //Some place to store the histograms
   std::vector<TH1*> PredictionHistograms;
@@ -42,64 +31,48 @@ int main(int argc, char * argv[]) {
   auto OutputFile = std::unique_ptr<TFile>(TFile::Open(OutputFileName.c_str(), "RECREATE"));
   OutputFile->cd();
 
-  osc->SetParameters(FitManager->raw()["General"]["OscillationParameters"].as<std::vector<double>>());
-  for (unsigned sample_i = 0 ; sample_i < DUNEPdfs.size() ; ++sample_i) {
-
-    std::string name = DUNEPdfs[sample_i]->GetTitle();
-    sample_names.push_back(name);
-    TString NameTString = TString(name.c_str());
+  for (auto handler : samples) {
+    for (unsigned iSample = 0; iSample < handler->GetNSamples(); ++iSample) {
 
-    DUNEPdfs[sample_i] -> Reweight();
-    if (DUNEPdfs[sample_i]->GetNDim() == 1){
-      PredictionHistograms.push_back(static_cast<TH1*>(DUNEPdfs[sample_i] -> Get1DHist() -> Clone(NameTString+"_unosc")));
-      DUNEPdfs[sample_i]->AddData(static_cast<TH1D*>(PredictionHistograms[sample_i]));
-    }
+      std::string name = handler->GetSampleTitle(iSample);
+      sample_names.push_back(name);
+      TString NameTString = TString(name.c_str());
 
-    else if (DUNEPdfs[sample_i]->GetNDim() == 2){
-      PredictionHistograms.push_back(static_cast<TH1*>(DUNEPdfs[sample_i] -> Get2DHist() -> Clone(NameTString+"_unosc")));
-      DUNEPdfs[sample_i]->AddData(static_cast<TH2D*>(PredictionHistograms[sample_i]));
-    }
-    else {
-      MACH3LOG_ERROR("Unsupported number of dimensions > 2 - Quitting"); 
-      throw MaCh3Exception(__FILE__ , __LINE__ );
+      handler->Reweight();
+      PredictionHistograms.push_back(static_cast<TH1*>(handler->GetMCHist(iSample)->Clone(NameTString+"_DataHist")));
+
+      if (handler->GetNDim(iSample) == 1){
+        handler->AddData(iSample, static_cast<TH1D*>(PredictionHistograms.back()));
+      } else if (handler->GetNDim(iSample) == 2){
+        handler->AddData(iSample, static_cast<TH2D*>(PredictionHistograms.back()));
+      }
+
+      else {
+        MACH3LOG_ERROR("Unsupported number of dimensions > 2 - Quitting"); 
+        throw MaCh3Exception(__FILE__ , __LINE__ );
+      }
+
+      MACH3LOG_INFO("Integrals of nominal hists: ");
+      MACH3LOG_INFO("{} : {}",name.c_str(),PredictionHistograms.back()->Integral());
+      MACH3LOG_INFO("--------------");
     }
-
-
-  }
-
-  //Now print out some event rates, we'll make a nice latex table at some point 
-  for (unsigned iPDF = 0; iPDF < DUNEPdfs.size() ; ++iPDF) {
-    MACH3LOG_INFO("Integrals of nominal hists: ");
-    MACH3LOG_INFO("{} : {}",sample_names[iPDF].c_str(),PredictionHistograms[iPDF]->Integral());
-    MACH3LOG_INFO("--------------");
   }
 
   //###########################################################################################################
   //MCMC
 
-  std::unique_ptr<mcmc> MaCh3Fitter = std::make_unique<mcmc>(FitManager);
+  auto MaCh3Fitter = MaCh3FitterFactory(FitManager.get());
 
   bool StartFromPreviousChain = GetFromManager(FitManager->raw()["General"]["StartFromPos"], false);
-
   //Start chain from random position unless continuing a chain
   if(!StartFromPreviousChain){
     if (!GetFromManager(FitManager->raw()["General"]["StatOnly"], false)) {
-      xsec->ThrowParameters();
+      param_handler->ThrowParameters();
     }
-    osc->SetParameters();
-    osc->ThrowParameters();
   }
-
 
   //Add systematic objects
-  MaCh3Fitter->AddSystObj(osc);
-  if (GetFromManager(FitManager->raw()["General"]["StatOnly"], false)){
-    MACH3LOG_INFO("Running a stat-only fit so no systematics will be applied");
-  }
-  else {
-    MaCh3Fitter->AddSystObj(xsec);
-  }
-
+  MaCh3Fitter->AddSystObj(param_handler.get());
 
   if (StartFromPreviousChain) {
     std::string PreviousChainPath = FitManager->raw()["General"]["PosFileName"].as<std::string>();
@@ -108,7 +81,7 @@ int main(int argc, char * argv[]) {
   }
 
   //Add samples
-  for(auto Sample : DUNEPdfs){
+  for(auto Sample : samples){
     MaCh3Fitter->AddSampleHandler(Sample);
   }
 
@@ -117,8 +90,8 @@ int main(int argc, char * argv[]) {
   MaCh3Fitter->RunMCMC();
 
   //Writing the memory usage at the end to eventually spot some nasty leak
-  MACH3LOG_WARN("\033[0;31mCurrent Total RAM usage is {:.2f} GB\033[0m", MaCh3Utils::getValue("VmRSS") / 1048576.0);
-  MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", MaCh3Utils::getValue("MemTotal") / 1048576.0);
+  MACH3LOG_WARN("\033[0;31mCurrent Total RAM usage is {:.2f} GB\033[0m", M3::Utils::getValue("VmRSS") / 1048576.0);
+  MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", M3::Utils::getValue("MemTotal") / 1048576.0);
 
   return 0;
 }