Problem implementing parallel fill of an array after ~4 hours of the code running (C++/ROOT)

Question

I am trying to implement parallel processing to fill an array with results from a function that runs an mcmc. Every test (using less loops in the chain, and less parameters to fit) I've done indicates that everything is working, but for some reason I get the following error after ~4 hours of the code running (this has happened twice):

root.exe(821,0x700001891000) malloc: *** error for object 0x6000015f3690: pointer being freed was not allocated root.exe(821,0x700001685000) malloc: *** error for object 0x6000015f3600: pointer being freed was not allocated

Here is the snippet where the parallel processing is being implemented (please note that filling the array with the mh_mcmc function worked perfectly before, so I don't think that is the issue). I can post more of the code as well if that would help:

Mat<double> mcmc_results[paramBins][paramBins][paramBins]; // array to store output of MCMC
// now we fill the mcmc_results array using multiple threads
// parallelization follows this: https://root.cern/doc/master/mt201__parallelHistoFill_8C.html`

gROOT->SetBatch(); // Avoid unnecessary output
ROOT::EnableThreadSafety(); // Make ROOT thread-safe
const int nWorkers = 4; // Match number of threads to what the hardware can do (my laptop has 4 cores

auto workFunction = [&](int k){
    
    arma::Mat<double> paramVals(paramBins, 3); paramVals = getParamVals();
    arma::Col<double> t23_col = paramVals.col(0);
    arma::Col<double> dmsq31_col = paramVals.col(1);
    arma::Col<double> dCP_col = paramVals.col(2);
    
    arma::Col<double> E_nu = getEnergyVals();
    
    for(int i=0; i<paramBins; i++){ // theta23
        cout << "i = " << i << endl;
        for(int j=0; j<paramBins; j++){ // dmsq31
            cout << "j = " << j << endl;
            
            arma::Col<double> params_obs = {t23_col(i),dmsq31_col(j),dCP_col(k)};
            mcmc_results[i][j][k] = mh_mcmc(params_obs, params_obs, E_nu, nu_flavor, NuOrAnti, muon_polarization, N_loops, param_cov);
            
        }
    }
};


// Create worker threads
std::vector<std::thread> workers;
for (auto k : ROOT::TSeqI(paramBins)){
    workers.emplace_back(workFunction, k);
}
// Make sure workers are done
for (auto &&worker : workers){
    worker.join();
}

I am a physicist, NOT a programmer (this will likely be the first and last time I ever have to do this sort of thing) so I am not sure what is going wrong. I have tried googling my problem, but haven't seen anything similar. Could this mean that my laptop isn't capable of storing all of this information (like not enough RAM or something)? Or could someone point me in the correct direction?

EDIT: The code worked perfectly fine before I tried to implement parallelization, so I dont think the problem is in the mh_mcmc function (the function performs a metropolis hastings mcmc). Here is the entire code:

    ''''
//
//  MH_MCMC.cpp
//  
//
//  Created by Grace Reesman on 12/5/22.
//

#include "NeutrinoFactory_MCMC.hpp"
#include "/Users/greesman/Desktop/NeutrinoFactory/Code/NeutrinoFactory.hpp"

#include <stdio.h>
#include <cmath> // math stuff
#include <complex>
#include <algorithm>
#include <cassert> // use the "assert" function
#include <cstddef>
#include <tgmath.h>
#include <vector> // to define std::vec
#include <iostream> // for cout/endl stuff (i think)
#include <fstream> // for cout/endl stuff (i think)
#include </opt/local/include/armadillo>
#include <limits> // this is to "define" infinity
//#include </Users/greesman/opt/anaconda3/pkgs/llvm-openmp-12.0.0-h0dcd299_1/include/omp.h> // parallel processing

// ROOT stuff to include
#include <TH3D.h> // 2D histograms
#include <TH2D.h> // 2D histograms
#include <TH1D.h> // 1D histograms
#include <TF1.h> // 1D functions
#include <TF2.h> // 2D functions
#include <TRandom3.h> // random numbers
#include <TLatex.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TFile.h>
#include <TGraph.h> // 1D graphs
#include <TGraph2D.h> // 2D graphs
#include <TFrame.h>
//#include <ROOT/TProcessExecutor.hxx> // parallel processing
#include <ROOT/TSeq.hxx> // for sequences
//#include <ROOT/TExecutor.hxx> // executing parallel processing
//#include <TThread.h>
#include <ROOT/RDataFrame.hxx>
#include <ROOT/RVec.hxx>


using namespace std;
using namespace arma;

Mat<double> getParamVals(){
    //std::cout << "Getting column vectors of parameter values..." << std::endl;
    // defining column vectors that will have paramBins entries each, where the values are uniform over the 3sigma nuFit estimates
    Col<double> dCP_col(paramBins, fill::zeros), t23_col(paramBins, fill::zeros), dmsq31_col(paramBins, fill::zeros);
    
    // setting minimum values
    dCP_col(0) = deltaCP_min;
    t23_col(0) = t23_min;
    dmsq31_col(0) = dmsq31_min;
    
    for(int i=1; i<paramBins; i++){ // filling the column vectors
        dCP_col(i) = dCP_col(i-1) + (deltaCP_max-deltaCP_min)/paramBins;
        t23_col(i) = t23_col(i-1) + (t23_max-t23_min)/paramBins;
        dmsq31_col(i) = dmsq31_col(i-1) + (dmsq31_max-dmsq31_min)/paramBins;
    }
    
    Mat<double> paramVals(paramBins, 3); // define a matrix that is MxN where M = number of parameters and N = number of param bins
    paramVals.col(0) = t23_col;
    paramVals.col(1) = dmsq31_col;
    paramVals.col(2) = dCP_col;
    
    return paramVals;
}

Col<double> getEnergyVals(){
    //std::cout << "Getting column vector of neutrino energies..." << std::endl;
    Col<double> E_nu(energyBins, fill::zeros); // defining and initializing a column vector to store energy values
    E_nu(0) = minEnergy; // setting the starting energy to the minimum energy
    for(int i = 1; i<energyBins; i++){
        E_nu(i) = E_nu(i-1) + step_size; // writing to neutrino energy column
    }
    return E_nu;
}

double radians(double y){ // degrees to radians
    return y * (M_PI / 180);
}

double normal_pdf(double x, double sigma, double x0) {
       double tmp = (x-x0)/sigma;
       return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
    }

// function that returns the model predictions
double get_model_predictions(Col<double> theta, double energy, int nu_flavor, int NuOrAntiNu, int mu_pol){
    // theta contains values of oscillation parameters of interest: theta(0) = theta23, theta(1) = delta m^2_{31}, theta(2) = delta_{CP}
    // energy is neutrino energy
    // nu_flavor is the neutrino flavor (or rather, the corresponding lepton) we want to look for in the far detector
    // NuOrAntiNu is 1 if neutrino, -1 if anti neutrino
    // mu_pol is muon polarization
    
    assert (nu_flavor >= 0 && nu_flavor <= 2); // nu_flavor valid for 0 (electron flavored), 1 (muon flavored), 2 (tau flavored)
    assert (NuOrAntiNu == -1 || NuOrAntiNu == 1); // NuOrAntiNu valid for -1 (anti-neutrino), 1 (neutrino)
    assert (mu_pol >= -1 && mu_pol <= 1); // muon polarization valid between -1 and 1
    
    double matter_potential = NuOrAntiNu*Yrho*Y_to_a*energy;
    
    // calculate the expected flux each neutrino flavor
    double ExpectFlux_NuE = Flux::diffFlux(0, NuOrAntiNu,mu_pol,energy)*Exact::Palphabeta(0,nu_flavor,matter_potential,L,energy,t12,theta(0),t13,Dmsq21,theta(1),NuOrAntiNu*theta(2));
    double ExpectFlux_NuMu = Flux::diffFlux(1, NuOrAntiNu, mu_pol, energy)*Exact::Palphabeta(1,nu_flavor,matter_potential,L,energy,t12,theta(0),t13,Dmsq21,theta(1),NuOrAntiNu*theta(2));
    double ExpectFlux_NuTau = Flux::diffFlux(2, NuOrAntiNu, mu_pol, energy)*Exact::Palphabeta(2,nu_flavor,matter_potential,L,energy,t12,theta(0),t13,Dmsq21,theta(1),NuOrAntiNu*theta(2)); // should be zero
    
    double ExpectCCXSec = CrossSection::sigma_nu(energy, NuOrAntiNu);
    return ExpectCCXSec*(ExpectFlux_NuE+ExpectFlux_NuMu+ExpectFlux_NuTau);
}

// function that returns the log of the likelihood
double lnlikelihood(Col<double> params_theory, Col<double> params_obs, Col<double> x_energy, int nu_flavor, int NuOrAntiNu, int mu_pol){
    
    double lnl = 0; // initializing the log likelihood
    
    assert(params_theory.n_elem == params_obs.n_elem); // make sure the two column vectors containing the parameters are the same size
    
    for(uword i=0; i<x_energy.n_elem; i++){
        
        double model_pred_theory = get_model_predictions(params_theory,x_energy(i),nu_flavor,NuOrAntiNu,mu_pol);
        double model_pred_obs = get_model_predictions(params_obs,x_energy(i),nu_flavor,NuOrAntiNu,mu_pol);
        
        lnl +=  -0.5*pow((model_pred_obs-model_pred_theory),2)/pow(sqrt(model_pred_obs),2); // add term to log likelihood
    }
    
    return lnl; // return the log of the likelihood
}

// function that calculates the log of the prior
double lnprior(Col<double> theta){
    
    double t23_input = theta(0);
    double dmsq31_input = theta(1);
    double dCP_input = theta(2);
    
    double lnp;
    
    double t23_prior = normal_pdf(t23_input,3*t23_1sig,t23);
    double dmsq31_prior = normal_pdf(dmsq31_input,3*dmsq31_1sig,Dmsq31);
    //+0.3*normal_pdf(-dmsq31_input,3*dmsq31_1sig,Dmsq31); // this is broken up into two parts so that we can see if we get any info about inverted ordering
    double dCP_prior = normal_pdf(dCP_input,3*dCP_1sig,deltaCP);
    
    if((t23_input<t23_zeroBelow)||(t23_input>t23_zeroAbove)||(dmsq31_input<dmsq31_zeroBelow)||(dmsq31_input>dmsq31_zeroAbove)||(dCP_input<dCP_zeroBelow)||(dCP_input>dCP_zeroAbove)){
        lnp = -infty;
    }
    else{
        lnp = log(t23_prior)+log(dmsq31_prior)+log(dCP_prior);
    }
    return lnp;
}

// test function that calculates the log of the unnormalized posterior
double lnposterior(Col<double> params_theory, Col<double> params_obs, Col<double> x_energy, int nu_flavor, int NuOrAntiNu, int mu_pol){
    
    double lnl = lnlikelihood(params_theory,params_obs,x_energy,nu_flavor,NuOrAntiNu,mu_pol);
    double lnp = lnprior(params_theory);
    double lnpost;
    
    if((lnp<=-infty/(10e15))||(lnl<=-infty/(10e15))){
        lnpost = -infty;
    }
    else{
        lnpost = lnl + lnp;
    }
    return lnpost;
}

// function to calculate the Hastings ratio
double hastings_ratio(Col<double> params_theory1, Col<double> params_theory2, Col<double> params_obs, Col<double> x_energy, int nu_flavor, int NuOrAntiNu, int mu_pol){
    
    double lnpost1 = lnposterior(params_theory1, params_obs, x_energy, nu_flavor, NuOrAntiNu, mu_pol);
    double lnpost2 = lnposterior(params_theory2, params_obs, x_energy, nu_flavor, NuOrAntiNu, mu_pol);
    double h_ratio;
    
    if((lnpost1 <= -infty/(10e15)) || (lnpost2 <= -infty/(10e15))){
        h_ratio =  -infty;
    }
    else{
        h_ratio = exp(lnpost1-lnpost2);
    }
    //cout << "ln posteriors to compare: " << lnpost1 << "and" << lnpost2 << endl;
    return h_ratio;
}


// function to generate a proposed new position for MCMC chain
Mat<double> propose_jump(Col<double> theta, Mat<double> cov){
    
    arma::arma_rng::set_seed_random();
    const int N = theta.n_elem;
    Col<double> proposed_position(N);
    proposed_position= mvnrnd(theta, cov); // generate a matrix w/ random column vectors from a multivariate Gaussian (normal) distribution
    
    while((proposed_position(0)<t23_zeroBelow)||(proposed_position(0)>t23_zeroAbove)||(proposed_position(1)<dmsq31_zeroBelow)||(proposed_position(1)>dmsq31_zeroAbove)||(proposed_position(2)<dCP_zeroBelow)||(proposed_position(2)>dCP_zeroAbove)){
        proposed_position = mvnrnd(theta, cov);
    }
    return proposed_position;
}

Mat<double> mh_mcmc(Col<double> params_theory, Col<double> params_obs, Col<double> x_energy, int nu_flavor, int NuOrAntiNu, int mu_pol, int N_steps, Mat<double> cov){
    //arma::arma_rng::set_seed_random();
    
    // here we will define/initialize some matrices/col vecs that we want to return
    Mat<double> positions(N_steps+1,params_theory.n_elem); // define a matrix that is MxN where M = number of steps to take in MCMC chain and N = number of model parameters
    Col<double> lnpost_at_pos(N_steps+1); // column vector that stores log-posterior value at the position of the MCMC chain
    Col<double> acceptance_ratio(N_steps+1); // column vector that stores acceptance ratio of all previous steps in the chain
    
    positions.row(0) = trans(params_theory);
    lnpost_at_pos(0) = lnposterior(params_theory,params_obs,x_energy,nu_flavor,NuOrAntiNu,mu_pol);
    int accepted = 0;
    
    for(uword step=0; step<N_steps; step++){
        //bool check = FALSE;
        
        Col<double> theta_propose = propose_jump(trans(positions.row(step)),cov);
        //cout << "proposed step: " << endl; theta_propose.print();
        
        double hastings = hastings_ratio(theta_propose, trans(positions.row(step)), params_theory, x_energy, nu_flavor,NuOrAntiNu,mu_pol);
        //cout << "Hastings ratio for proposed step: " << hastings << endl;

        Col<double> rand_colvec(1,fill::randu); double p_rand = rand_colvec(0); // generating a random number in (0,1)
        //cout << "Random # to compare to hastings: " << p_rand << endl;

        if(hastings>=p_rand){
            positions.row(step+1) = trans(propose_jump(trans(positions.row(step)),cov));
            lnpost_at_pos(step+1) = lnposterior(trans(positions.row(step+1)),params_obs,x_energy,nu_flavor,NuOrAntiNu,mu_pol);
            accepted += 1;
            acceptance_ratio(step+1) = accepted/(step+1);
            //check = TRUE;
        }
        else{
            positions.row(step+1) = positions.row(step);
            lnpost_at_pos.row(step+1) = lnpost_at_pos(step);
            acceptance_ratio(step+1) = accepted/(step+1);
        }
        //cout << check << endl;
    }
    
    //cout << "Number of accepted positions: " << accepted << endl;
    Mat<double> returnMat(N_steps+1,5);
    returnMat.col(0) = positions.col(0);
    returnMat.col(1) = positions.col(1);
    returnMat.col(2) = positions.col(2);
    returnMat.col(3) = lnpost_at_pos;
    returnMat.col(4) = acceptance_ratio;
        
    return returnMat;
}

TGraph2D* mcmcPositions_graph(Mat<double> MCMC_results, int Nloops, int name1, int name2, int name3){
    assert(MCMC_results.n_rows == Nloops+1);
    
    Mat<double> MCMC_positions = MCMC_results.head_cols(3);
    
    double t23_mcmc_arr[Nloops+1], dmsq31_mcmc_arr[Nloops+1], dCP_mcmc_arr[Nloops+1];
    for(int i=0; i<Nloops+1; i++){
        t23_mcmc_arr[i] = 0;
        dmsq31_mcmc_arr[i] = 0;
        dCP_mcmc_arr[i] = 0;
    }
    for(int i=0; i<Nloops+1; i++){
        t23_mcmc_arr[i] = MCMC_positions(i,0);
        dmsq31_mcmc_arr[i] = MCMC_positions(i,1);
        dCP_mcmc_arr[i] = MCMC_positions(i,2);
    }
    
    auto param_mcmc_gr = new TGraph2D(Nloops+1,t23_mcmc_arr, dmsq31_mcmc_arr, dCP_mcmc_arr);
    param_mcmc_gr->SetName(Form("MCMCPositions%d%d%d",name1,name2,name3));
    param_mcmc_gr->SetTitle(Form("MCMCPositions%d%d%d",name1,name2,name3));
    param_mcmc_gr->GetXaxis()->SetTitle("#theta_{23}");
    param_mcmc_gr->GetYaxis()->SetTitle("#Delta m^2_{31}");
    param_mcmc_gr->GetZaxis()->SetTitle("#delta_{CP}");
    param_mcmc_gr->SetMarkerColorAlpha(4,0.25);
    param_mcmc_gr->SetMarkerStyle(21);
    
    return param_mcmc_gr;
}

TH3D* mcmcPositions_hist(Mat<double> MCMC_results, int Nloops){
    
    TGraph2D* mcmc_gr = mcmcPositions_graph(MCMC_results,Nloops,0,0,0);
    TH3D* mcmc_hist = new TH3D("mcmc_hist", Form("%s;%s;%s;%s", mcmc_gr->GetTitle(), mcmc_gr->GetXaxis()->GetTitle(), mcmc_gr->GetYaxis()->GetTitle(), mcmc_gr->GetZaxis()->GetTitle()),
                               paramBins, t23_min, t23_max,
                               paramBins, dmsq31_min, dmsq31_max,
                               paramBins, deltaCP_min, deltaCP_max);
    
    for (int i=0; i < mcmc_gr->GetN(); i++){ // setting bin contents to the TGraph values
        double t23_pt,dmsq31_pt,dCP_pt;
        mcmc_gr->GetPoint(i,t23_pt,dmsq31_pt,dCP_pt);
        mcmc_hist->Fill(t23_pt,dmsq31_pt,dCP_pt);
    }
    mcmc_hist->Scale(1/mcmc_hist->Integral("width")); // normalizing our histogram
    return mcmc_hist;
}

TH2D* TwoParam_pdf(Mat<double> MCMC_results, int Nloops, int margeParam){
    
    assert(0<=margeParam && margeParam<=2); // margParam is the parameter we want to marginalize over (0 is t23, 1 is dmsq31, 2 is deltaCP)
    
    TH3D* ThreeParHist = mcmcPositions_hist(MCMC_results,Nloops);
    TH2D* TwoParamHist = new TH2D();
    
    if(margeParam==0){ //marginalize over t23
        TwoParamHist = static_cast<TH2D*>(ThreeParHist->Project3D("yz"));
    }
    else if(margeParam==1){ //marginalize over dmsq31
        TwoParamHist = static_cast<TH2D*>(ThreeParHist->Project3D("xz"));
    }
    else if(margeParam==2){ //marginalize over dCP
        TwoParamHist = static_cast<TH2D*>(ThreeParHist->Project3D("xy"));
    }
    else{
        cout << "problem marginalizing over one parameter" << endl;
    }
    return TwoParamHist;
}

TH1D* OneParam_pdf(Mat<double> MCMC_results, int Nloops, int param){
    
    assert(0<=param && param<=2); // param is the only parameter we won't marginalize over (0 is t23, 1 is dmsq31, 2 is deltaCP)
    
    TH3D* ThreeParHist = mcmcPositions_hist(MCMC_results,Nloops);
    TH1D* OneParamHist = new TH1D();
    
    if(param==0){ //marginalize over t23
        OneParamHist = static_cast<TH1D*>(ThreeParHist->Project3D("x"));
    }
    else if(param==1){ //marginalize over dmsq31
        OneParamHist = static_cast<TH1D*>(ThreeParHist->Project3D("y"));
    }
    else if(param==2){ //marginalize over dCP
        OneParamHist = static_cast<TH1D*>(ThreeParHist->Project3D("z"));
    }
    else{
        cout << "problem marginalizing over two parameters" << endl;
    }
    return OneParamHist;
}

void plot_mh_summary(Mat<double> MCMC_results, int Nloops){
    
    // one dimensional histograms
    TH1D* t23_hist = OneParam_pdf(MCMC_results,Nloops,0);
    TH1D* dmsq31_hist = OneParam_pdf(MCMC_results,Nloops,1);
    TH1D* dCP_hist = OneParam_pdf(MCMC_results,Nloops,2);
    
    // canvas
    auto OneD_canvas = new TCanvas("OneD_canvas","",200,10,700,500);
    //OneD_canvas->SetFillColor(42);
    OneD_canvas->SetGrid();
    OneD_canvas->GetFrame()->SetFillColor(21);
    OneD_canvas->GetFrame()->SetBorderSize(12);
    OneD_canvas->Divide(3);
    
    OneD_canvas->cd(1);
    t23_hist->Draw("hist");
    OneD_canvas->cd(2);
    dmsq31_hist->Draw("hist");
    OneD_canvas->cd(3);
    dCP_hist->Draw("hist");
    
    // two dimensional histograms
    TH2D* t23_dmsq31_hist = TwoParam_pdf(MCMC_results,Nloops,2);
    TH2D* t23_dCP_hist = TwoParam_pdf(MCMC_results,Nloops,1);
    TH2D* dmsq31_dCP_hist = TwoParam_pdf(MCMC_results,Nloops,0);
                                         
    // canvas
    auto TwoD_canvas = new TCanvas("TwoD_canvas","",200,10,700,500);
    //TwoD_canvas->SetFillColor(42);
    TwoD_canvas->SetGrid();
    TwoD_canvas->GetFrame()->SetFillColor(21);
    TwoD_canvas->GetFrame()->SetBorderSize(12);
    TwoD_canvas->Divide(3);
    
    TwoD_canvas->cd(1);
    t23_dmsq31_hist->Draw();
    TwoD_canvas->cd(2);
    t23_dCP_hist->Draw();
    TwoD_canvas->cd(3);
    dmsq31_dCP_hist->Draw();
    
}

TH1D* Col2TH1D(Col<double> col_vec, int name1, int name2, int name3){
    
    int N = col_vec.n_elem;
    double min_val = min(col_vec);
    double max_val = max(col_vec);
    TH1D* hist = new TH1D("hist","",N,min_val,max_val);
    hist->SetName(Form("MCMCPositions%d%d%d",name1,name2,name3));
    
    for(uword i=0; i<N; i++){
        hist->Fill(col_vec(i));
    }
    return hist;
}

Mat<double> mcmc_results[paramBins][paramBins][paramBins]; // array to store output of MCMC

// NOTE: this function assumes a muon beam (=> muon polarization = +1 and wrong sign muon signal will be antimuons)
void eventSpecific_mcmc(int lep_flavor, int lep_antilep){
    // variables indicate whether the lepton flavor and whether we have a lepton or antilepton
    // lep_flavor key: 0~electron, 1~muon, 2~tauon
    // lep_antilep key: 1~lepton, -1~anti-lepton
    
    assert (lep_flavor >= 0 && lep_flavor <= 2); // make sure lep_flavor is in the correct range
    assert (lep_antilep == 1 || lep_antilep == -1); // make sure lep_flavor is in the correct range
    
    string eventType; // variable to store the type of event we are doing the MCMC for in string form
    if(lep_flavor == 1){
        if(lep_antilep == -1){
            eventType = "wMu";
        }
        else if(lep_antilep == 1){
            eventType = "rMu";
        }
    }
    else if(lep_flavor == 0){
        if(lep_antilep == -1){
            eventType = "rElec";
        }
        else if(lep_antilep == 1){
            eventType = "wElec";
        }
    }
    else{
        if(lep_antilep == -1){
            eventType = "antiTau";
        }
        else if(lep_antilep == 1){
            eventType = "Tau";
        }
    }
    
    Mat<double> mcmc_results[paramBins][paramBins][paramBins]; // array to store output of MCMC
    // now we fill the mcmc_results array using multiple threads
    // parallelization follows this: https://root.cern/doc/master/mt201__parallelHistoFill_8C.html
    
    gROOT->SetBatch(); // Avoid unnecessary output
    ROOT::EnableThreadSafety(); // Make ROOT thread-safe
    const int nWorkers = 4; // Match number of threads to what the hardware can do (my laptop has 4 cores
    
    auto workFunction = [&](int k){
        
        Mat<double> paramVals(paramBins, 3); paramVals = getParamVals();
        Col<double> t23_col = paramVals.col(0);
        Col<double> dmsq31_col = paramVals.col(1);
        Col<double> dCP_col = paramVals.col(2);
        
        Col<double> E_nu = getEnergyVals();
        
        for(int i=0; i<paramBins; i++){ // theta23
            cout << "i = " << i << endl;
            for(int j=0; j<paramBins; j++){ // dmsq31
                cout << "j = " << j << endl;
                
                Col<double> params_obs = {t23_col(i),dmsq31_col(j),dCP_col(k)};
                mcmc_results[i][j][k] = mh_mcmc(params_obs, params_obs, E_nu, nu_flavor, NuOrAnti, muon_polarization, N_loops, param_cov);
                
            }
        }
    };
    
    
    // Create worker threads
    std::vector<std::thread> workers;
    for (auto k : ROOT::TSeqI(paramBins)){
        workers.emplace_back(workFunction, k);
    }
    // Make sure workers are done
    for (auto &&worker : workers){
        worker.join();
    }
    
    // files to save the output of MCMC
    TFile* Positions_outFile = new TFile(Form("/Users/greesman/Desktop/NeutrinoFactory/OutputFiles/Positions/%s/%s_mcmcPositions.root",eventType.c_str(),eventType.c_str()),"RECREATE");
    for(int i=0; i<paramBins; i++){
        for(int j=0; j<paramBins; j++){
            for(int k=0; k<paramBins; k++){
                mcmcPositions_graph(mcmc_results[i][j][k],N_loops,i,j,k)->Write();
            }
        }
    }

    TFile* Posteriors_outFile = new TFile(Form("/Users/greesman/Desktop/NeutrinoFactory/OutputFiles/Posteriors/%s/%s_mcmcPosteriors.root",eventType.c_str(),eventType.c_str()),"RECREATE");
    for(int i=0; i<paramBins; i++){
        for(int j=0; j<paramBins; j++){
            for(int k=0; k<paramBins; k++){
                Col2TH1D(mcmc_results[i][j][k].col(3),i,j,k)->Write();
            }
        }
    }

    TFile* AccRatio_outFile = new TFile(Form("/Users/greesman/Desktop/NeutrinoFactory/OutputFiles/AcceptanceRatios/%s/%s_mcmcAcceptanceRatios.root",eventType.c_str(),eventType.c_str()),"RECREATE");
    for(int i=0; i<paramBins; i++){
        for(int j=0; j<paramBins; j++){
            for(int k=0; k<paramBins; k++){
                Col2TH1D(mcmc_results[i][j][k].col(4),i,j,k)->Write();
            }
        }
    }

}

void NeutrinoFactory_MCMC(){
    
    std::cout << "Starting..." << std::endl;
    
    std::cout << "performing MH MCMC for wrong-sign muons" << std::endl;
    eventSpecific_mcmc(1,-1); //wMus
    //std::cout << "performing MH MCMC for right-sign muons" << std::endl;
    //eventSpecific_mcmc(1,1); //rMus
    //std::cout << "performing MH MCMC for wrong-sign electrons" << std::endl;
    //eventSpecific_mcmc(0,1); //wElec
    //std::cout << "performing MH MCMC for right-sign electrons" << std::endl;
    //eventSpecific_mcmc(0,-1); //rElec
    

}
''''