Merge pull request #43 from DUNE/atm_chains

Adding elements to use chains with atmospherics

configs/EventRates_Atmospherics.yaml

@@ -37,6 +37,8 @@ General:
   ProcessMCMC: No
   Seed: 0
   Debug: No
+  StartFromPos: false
+  PosFileName: "chain_0.root"
 
 "Projections": [
   {

make_tar.sh

@@ -0,0 +1 @@
+ tar --exclude='**/.git' -C .. -cvf ../mach3.tar.gz MaCh3_DUNE

setup_dune_env.sh

@@ -1,4 +1,4 @@
 source /cvmfs/larsoft.opensciencegrid.org/spack-packages/setup-env.sh
-spack load [email protected].06
+spack load [email protected].12
 spack load cmake
 spack load [email protected]

setup_grid.sh

@@ -0,0 +1,4 @@
+source ./build/bin/setup.MaCh3.sh
+source ./build/bin/setup.MaCh3DUNE.sh
+source ./build/bin/setup.NuOscillator.sh
+spack load [email protected]

src/Fit.cpp

@@ -36,7 +36,7 @@ int main(int argc, char * argv[]) {
   MakeMaCh3DuneInstance(FitManager, DUNEPdfs, xsec, osc); 
 
   //Some place to store the histograms
-  std::vector<TH1D*> PredictionHistograms;
+  std::vector<TH1*> PredictionHistograms;
   std::vector<std::string> sample_names;
 
   auto OutputFile = std::unique_ptr<TFile>(TFile::Open(OutputFileName.c_str(), "RECREATE"));
@@ -50,10 +50,21 @@ int main(int argc, char * argv[]) {
 
     osc->setParameters();
     DUNEPdfs[sample_i] -> reweight();
-    TH1D *SampleHistogram = (TH1D*)DUNEPdfs[sample_i] -> get1DHist() -> Clone(NameTString+"_unosc");
-    PredictionHistograms.push_back(SampleHistogram);
+    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]));
+    }
+
+    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__ );
+    }
+
 
-    DUNEPdfs[sample_i]->addData(PredictionHistograms[sample_i]);
   }
 
   //Now print out some event rates, we'll make a nice latex table at some point 
@@ -69,20 +80,13 @@ int main(int argc, char * argv[]) {
   std::unique_ptr<mcmc> MaCh3Fitter = std::make_unique<mcmc>(FitManager);
 
   bool StartFromPreviousChain = GetFromManager(FitManager->raw()["General"]["StartFromPos"], false);
-  if (StartFromPreviousChain) {
-    std::string PreviousChainPath = FitManager->raw()["General"]["PosFileName"].as<std::string>();
-    MACH3LOG_INFO("MCMC getting starting position from: {}",PreviousChainPath);
-    MaCh3Fitter->StartFromPreviousFit(PreviousChainPath);
+
+  //Start chain from random position unless continuing a chain
+  if(!StartFromPreviousChain){
+    xsec->throwParameters();
+    osc->throwParameters();
   }
 
-  //Add samples
-  for(auto Sample : DUNEPdfs){
-    MaCh3Fitter->addSamplePDF(Sample);
-  }
-
-  //Start chain from random position
-  xsec->throwParameters();
-  osc->throwParameters();
 
   //Add systematic objects
   MaCh3Fitter->addSystObj(osc);
@@ -92,9 +96,26 @@ int main(int argc, char * argv[]) {
   else {
     MaCh3Fitter->addSystObj(xsec);
   }
+
+
+  if (StartFromPreviousChain) {
+    std::string PreviousChainPath = FitManager->raw()["General"]["PosFileName"].as<std::string>();
+    MACH3LOG_INFO("MCMC getting starting position from: {}",PreviousChainPath);
+    MaCh3Fitter->StartFromPreviousFit(PreviousChainPath);
+  }
+
+  //Add samples
+  for(auto Sample : DUNEPdfs){
+    MaCh3Fitter->addSamplePDF(Sample);
+  }
+
 
   //Run fit
   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);
+
   return 0;
 }