/*
DB: 19/01/2021
----------------------
PrintEventRate.cpp

- Prints Atmsopheric event rates to console and to a text file
- Currently includes AtmFlux and Xsec
- Text file hardcoded in text
- Text file is in latex format

Usage:
./PrintEvenRate config

Output:
Event rates formatetted in Latex to hard-coded file name (Currently EventRate_UnOsc.txt)
----------------------
*/

#include <iostream>
#include <vector>
#include <sys/time.h>

#include <TH1D.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TStopwatch.h>

#include "samplePDF/samplePDFSKBase.h"
#include "samplePDF/samplePDFSKBaseGPU.h"
#include "samplePDF/samplePDFND2019.h"
#include "samplePDF/samplePDFND2019GPU.h"
#include "covariance/covariancePsyche.h"
#include "mcmc/mcmc.h"
#include "manager/manager.h"
#include "samplePDF/Structs.h"
#include "covariance/BANFF_to_xsec_mapping.h"


#if defined(CUDA)
void PrintIntegral(std::vector<samplePDFSKBaseGPU*> PDFS, int WeightStyle, bool IncMigration,TString OutputName, TString Iden);
#else
void PrintIntegral(std::vector<samplePDFSKBase*> PDFS, int WeightStyle, bool IncMigration,TString OutputName, TString Iden);
#endif

void PrintToScreen(TString SampleName,std::vector< std::vector<double> > IntegralList, std::vector<TString> OscChannelNames, std::vector<TString> ModeNames,TString Iden,TString OutputFileName="/dev/null");

int main(int argc, char **argv) {
  TApplication a("a", 0, 0);

  gStyle->SetOptStat(0);
  int WeightStyle = 0; // 0 Use Event Weight, 1 Count Events

  bool useT2K = false;
  bool IncMigration = true;

  //####################################################################################
  //Get FitMan parameters
  manager *fitMan = new manager(argv[1]);

  std::cout << "\n\n" << "********* Parameters pulled from config file in Exec *********" << "\n" << std::endl;

  double NuPOT = fitMan->getPOT();
  std::cout << "NuPOT: " << NuPOT << std::endl;
  double NuBarPOT = fitMan->getNubarPOT();
  std::cout << "NuBarPOT: " << NuBarPOT << std::endl;

  double AtmLivetime = fitMan->getLivetime();
  std::cout << "Atmospheric MC livetime (Years): " << AtmLivetime << std::endl;


  std::vector<int> AtmSamples = fitMan->getAtmSamples();
  std::cout << "**** Samples:";
  for (int i=0;i<AtmSamples.size();i++) {
    std::cout << AtmSamples[i] << " ";
  }
  std::cout << " ****" << std::endl;

  std::vector<double> oscpars = fitMan->getOscParameters();
  std::cout << "Oscillation Parameters: { ";
  for (unsigned int i=0;i<oscpars.size();i++) {
    std::cout << oscpars[i] << " ";
  }
  std::cout << "}" << std::endl;

  TString oscCovMatrixFile  = fitMan -> getOscCovMatrix();
  TString oscCovMatrixName  = fitMan -> getOscCovName();
  std::cout << "Oscillation covariance used: " << oscCovMatrixName << " from " << oscCovMatrixFile << std::endl;

  std::string AtmFluxSystCfgName = fitMan->GetAtmFluxSystCfgName();
  std::cout << "Config used for Atmospheric Flux Systematics: " << AtmFluxSystCfgName << std::endl;

  std::string OscillatorCfgName = fitMan->GetOscillatorCfgName();
  std::cout << "Config used for Oscillator: " << OscillatorCfgName << std::endl;

  std::string SKCalibrationSystCfgName = fitMan->GetSKCalibrationSystCfgName();
  std::cout << "Config used for SK Calibration Systematics: " << SKCalibrationSystCfgName << std::endl;

  std::string SKDetBeamSystCfgName = fitMan->GetSKDetBeamSystCfgName();
  std::cout << "Config used for SK Detector Beam Systematics: " << SKDetBeamSystCfgName << std::endl;

  std::string AtmSKDetSystCfgName = fitMan->GetAtmSKDetSystCfgName();
  std::cout << "Config used for Atmospheric Flux Systematics: " << AtmSKDetSystCfgName << std::endl;

  std::string ATMPDDetectorSystCfgName = fitMan->GetATMPDDetectorSystCfgName();
  std::cout << "Config used for ATMPD Detector Systematics: " << ATMPDDetectorSystCfgName << std::endl;

  TString skDetCovMatrixFile = fitMan -> getSKdetCovMatrix();
  TString skDetCovMatrixName = fitMan -> getSKdetCovName();
  std::cout << "skDet covariance used: " << skDetCovMatrixName << " from " << skDetCovMatrixFile << std::endl;

  bool useBANFFTune = fitMan->getUseBANFFtuning();
  bool useGeneratedTune = fitMan->getUseGeneratedtuning();

  TString xsecCovMatrixFile  = fitMan -> getXsecCovMatrix();
  TString xsecCovMatrixName  = fitMan -> getXsecCovName();

  TString BANFFCovMatrixFile = fitMan -> getBANFFCovMatrix();
  TString BANFFCovMatrixName = fitMan -> getBANFFCovName() ;
  TString BANFFParamName     = fitMan -> getBANFFParamName();

  TString GeneratedTuneFile      = fitMan -> getGenTuneFile();
  TString GeneratedTuneParamName = fitMan -> getGenTuneParamName();

  std::cout << "Xsec covariance used: " << xsecCovMatrixName << " from " << xsecCovMatrixFile << std::endl;
  if (useBANFFTune && useGeneratedTune) {
    std::cout << "useBANFFTune and useGeneratedTune both true. Quitting.." << std::endl;
    throw;
  }
  if (useBANFFTune && !useGeneratedTune) {
    std::cout << "BANFF tune used: " << BANFFCovMatrixName << " from " << BANFFCovMatrixFile << std::endl;
  }
  if (!useBANFFTune && useGeneratedTune) {
    std::cout << "Generated tune used: " << GeneratedTuneParamName << " from " << GeneratedTuneFile << std::endl;
  }

  double xsecCovStepScale = fitMan->getXsecStepScale();
  std::cout << "Xsec covariance step scale: " << xsecCovStepScale << std::endl;

  double oscCovStepScale = fitMan->getOscStepScale();
  std::cout << "Osc covariance step scale: " << oscCovStepScale << std::endl;

  double T2KSKDetStepScale = fitMan->getSkDetStepScale();
  std::cout << "Beam SKDet covariance step scale: " << T2KSKDetStepScale << std::endl;

  double SKDetBeamStepScale = fitMan->getSKDetBeamStepScale();
  std::cout << "Additional SKDetBeam covariance step scale: " << SKDetBeamStepScale << std::endl;

  double AtmSKDetStepScale = fitMan->getAtmDetStepScale();
  std::cout << "AtmSKDet covariance step scale: " << AtmSKDetStepScale << std::endl;

  double ATMPDDetectorStepScale = fitMan->getATMPDDetStepScale();
  std::cout << "ATMPD Detector covariance step scale: " << ATMPDDetectorStepScale << std::endl;

  double SKCalibrationStepScale = fitMan->getSKCalibrationStepScale();
  std::cout << "SK Calibration covariance step scale: " << SKCalibrationStepScale << std::endl;

  double AtmFluxStepScale = fitMan->getAtmFluxStepScale();
  std::cout << "Atm Flux covariance step scale: " << AtmFluxStepScale << std::endl;

  std::vector<int> XsecFlatParams  = fitMan->getXsecFlat();
  std::cout << "Flat Xsec parameters: { ";
  for (unsigned int i=0;i<XsecFlatParams.size();i++) {
    std::cout << XsecFlatParams[i] << " ";
  }
  std::cout << "}" << std::endl;

  bool useRC = fitMan->getRC();
  if (useRC) {
    std::cout << "Using Reactor Constraint" << std::endl;
  } else {
    std::cout << "Not Using Reactor Constraint" << std::endl;
  }

  std::cout << "\n\n" << "******************" << "\n" << std::endl;

  //####################################################################################
  //Create Oscillator

  std::cout << "-----------------------------------------------------------" << std::endl;
  std::cout << "Creating Oscillator object" << std::endl;

  Oscillator* Oscill = new Oscillator(OscillatorCfgName);

  std::cout << "-----------------------------------------------------------" << std::endl;

  //####################################################################################
  //Set Osc

  std::vector<double> oscpars_un(oscpars);
  for (int i=0;i<3;i++) {
    oscpars_un[i] = 0.;
  }

  //Hardcoded 2020
  std::vector<double> oscpars_AsimovA(oscpars);
  oscpars_AsimovA [0] = 3.07e-1;
  oscpars_AsimovA [1] = 5.28e-1;
  oscpars_AsimovA [2] = 2.18e-2;
  oscpars_AsimovA [3] = 7.53e-5;
  oscpars_AsimovA [4] = 2.509e-3;
  oscpars_AsimovA [5] = -1.601;

  //Hardcoded 2020
  std::vector<double> oscpars_AsimovB(oscpars);
  oscpars_AsimovB [0] = 3.07e-1;
  oscpars_AsimovB [1] = 4.5e-1;
  oscpars_AsimovB [2] = 2.18e-2;
  oscpars_AsimovB [3] = 8.0e-5;
  oscpars_AsimovB [4] = 2.51e-3;
  oscpars_AsimovB [5] = 0.;

  covarianceOsc *osc = new covarianceOsc(oscCovMatrixName,oscCovMatrixFile);

  // Use prior for 12 parameters only
  //osc->setEvalLikelihood(0,false);
  osc->setEvalLikelihood(1,false);
  osc->setEvalLikelihood(2,false);
  //osc->setEvalLikelihood(3,false);
  osc->setEvalLikelihood(4,false);
  osc->setEvalLikelihood(5,false);

  osc->setFlipDeltaM23(true);
  osc->useReactorPrior(useRC);

  osc->setStepScale(oscCovStepScale);
  osc->setParameters(oscpars);
  osc->acceptStep();

  //####################################################################################
  //Set Atmospheric Flux

  covarianceAtmFlux* AtmFlux = new covarianceAtmFlux(AtmFluxSystCfgName);
  AtmFlux->setStepScale(AtmFluxStepScale);

  AtmFlux->setParameters();
  AtmFlux->printNominalCurrProp();

  //####################################################################################
  //Set Atmospheric SK Calibration Systematics

  covarianceSKCalibration* SKCalibration = new covarianceSKCalibration(SKCalibrationSystCfgName);
  SKCalibration->SetNominalEventRates();

  SKCalibration->setStepScale(SKCalibrationStepScale);

  SKCalibration->setParameters();
  SKCalibration->printNominalCurrProp();

  //####################################################################################
  //Set T2K OA SK Detector Systematics

  covarianceSkDet_joint* SKDet = new covarianceSkDet_joint(skDetCovMatrixName, skDetCovMatrixFile);
  SKDet->setStepScale(T2KSKDetStepScale);

  SKDet->setParameters();
  SKDet->printNominalCurrProp();

  //####################################################################################
  //Set T2K Supplementary Detector Systematics

  covarianceSKDetBeam* SKDetBeam = new covarianceSKDetBeam(SKDetBeamSystCfgName);
  SKDetBeam->setStepScale(SKDetBeamStepScale);

  SKDetBeam->setParameters();
  SKDetBeam->printNominalCurrProp();

  //####################################################################################
  //Set Atmospheric Detector

  covarianceAtmSKDet* AtmDet = new covarianceAtmSKDet(AtmSKDetSystCfgName);
  AtmDet->DoMigration(IncMigration);

  AtmDet->setStepScale(AtmSKDetStepScale);

  AtmDet->setParameters();
  AtmDet->printNominalCurrProp();

  //####################################################################################
  //Set ATMPD Detector

  covarianceATMPDDetector* ATMPDDet = new covarianceATMPDDetector(ATMPDDetectorSystCfgName);
  ATMPDDet->SetNominalEventRates();

  ATMPDDet->setStepScale(ATMPDDetectorStepScale);

  ATMPDDet->setParameters();
  ATMPDDet->printNominalCurrProp();

  //####################################################################################
  //Set Xsec
  covarianceXsec* xsec = new covarianceXsec(xsecCovMatrixName, xsecCovMatrixFile);

  if (XsecFlatParams.size() == 1 && XsecFlatParams.at(0) == -1) {
    for (int j = 0; j < xsec->getSize(); j++) {
      xsec->setEvalLikelihood(j, false);
    }
  } else {
    for (int j = 0; j < XsecFlatParams.size(); j++) {
      xsec->setEvalLikelihood(XsecFlatParams.at(j), false);
    }
  }

  xsec->setStepScale(xsecCovStepScale);

  //######################################
  // Find Pre and Post BANFF tunes

  std::vector<double> xsec_preBANFF = xsec->getNominalArray();
  std::vector<double> xsec_postBANFF = xsec->getNominalArray();
  std::vector<double> xsec_generated = xsec->getNominalArray();

  if (useBANFFTune) {
    //Get postfit parameter values from BANFF matrix
    TFile *BANFF_File = new TFile(BANFFCovMatrixFile,"READ");
    TVectorD *BANFF_ParamArray = (TVectorD*)BANFF_File->Get(BANFFParamName);

    BANFFtoMaCh3XsecCov_2020_Atm(xsec_postBANFF,*BANFF_ParamArray);

    std::cout << std::setw(35) << "Parameter:  " << " | " << std::setw(20) << "Pre-BANFF Value" << " | " << std::setw(20) << "Post-BANFF value" << std::endl;
    std::cout << "-------------------------------------------------------------------------------------" << std::endl;
    for (int iParam=0;iParam<xsec_preBANFF.size();iParam++) {
      std::cout << std::setw(35) << xsec->getParName(iParam) << " | " << std::setw(20) << xsec_preBANFF[iParam] << " | " << std::setw(20) << xsec_postBANFF[iParam] << std::endl;
    }
  }

  if (useGeneratedTune) {
    TFile* Xsec_File = new TFile(GeneratedTuneFile,"READ");
    TVectorD* Xsec_GeneratedParamArray = (TVectorD*)Xsec_File->Get(GeneratedTuneParamName);

    for(int iParam=0;iParam<xsec_generated.size();iParam++) {
      xsec_generated[iParam] = (*Xsec_GeneratedParamArray)(iParam);
    }
    std::cout << std::setw(35) << "Parameter:  " << " | " << std::setw(20) << "Pre-BANFF Value" << " | " << std::setw(20) << "Generated value" << std::endl;
    std::cout << "-------------------------------------------------------------------------------------" << std::endl;
    for (int iParam=0;iParam<xsec_preBANFF.size();iParam++) {
      std::cout << std::setw(35) << xsec->getParName(iParam) << " | " << std::setw(20) << xsec_preBANFF[iParam] << " | " << std::setw(20) << xsec_generated[iParam] << std::endl;
    }
  }

  TString Iden;

  if (useBANFFTune) {
    std::cout << "Setting Post BANFF Tune.." << std::endl;
    xsec->setParameters(xsec_postBANFF);
    Iden = "PostBANFF";
  } else if (useGeneratedTune) {
    std::cout << "Setting MC Generated Tune.." << std::endl;
    xsec->setParameters(xsec_generated);
    Iden = "Generated";
  } else {
    std::cout << "Setting Pre-BANFF Tune.." << std::endl;
    xsec->setParameters(xsec_preBANFF);

    xsec->setParProp(10,1.);
    xsec->setParProp(11,1.);
    xsec->setParProp(12,1.);
    xsec->setParProp(13,1.);
    xsec->setParProp(14,1.);

    Iden = "PreBANFF";
  }

  xsec->acceptStep();
  xsec->printNominalCurrProp();

  //####################################################################################
  //Create samplePDFSKBase Objs

#if defined(CUDA)
  std::vector<samplePDFSKBaseGPU*> AtmPdfs;
#else
  std::vector<samplePDFSKBase*> AtmPdfs;
#endif

  int nAtmSamplesLoaded= 0;

  for (int i=0;i<AtmSamples.size();i++) {

    TString cfgName = Form("configs/AtmosphericConfigs/AtmSample_%i.cfg",AtmSamples[i]);

#if defined(CUDA)
    samplePDFSKBaseGPU* ATM_PDF =  new samplePDFSKBaseGPU(AtmLivetime,cfgName.Data(),xsec);
#else
    samplePDFSKBase* ATM_PDF =  new samplePDFSKBase(AtmLivetime,cfgName.Data(),xsec);
#endif

    ATM_PDF->SetOscillator(Oscill);
    ATM_PDF->setUseBeta(false);
    ATM_PDF->setApplyBetaNue(false);
    ATM_PDF->setApplyBetaDiag(false);
    ATM_PDF->useNonDoubledAngles(true);
    std::cout << "-----------------------------" << std::endl;
    std::cout << "Setting Xsec covariance object.." << std::endl;
    ATM_PDF->setSkXsecCov(xsec);
    std::cout << "-----------------------------" << std::endl;
    std::cout << "Setting AtmFlux covariance object.." << std::endl;
    ATM_PDF->setAtmFluxCov(AtmFlux);
    std::cout << "-----------------------------" << std::endl;
    std::cout << "Setting SK Calibration covariance object.." << std::endl;
    ATM_PDF->setSKCalibrationCov(SKCalibration);

    if (ATM_PDF->GetSKDetObj() == kDetObj_AtmSKDet) {
      std::cout << "-----------------------------" << std::endl;
      std::cout << "Setting AtmSKDet covariance object.." << std::endl;
      ATM_PDF->setAtmSKDetCov(AtmDet);
    } else if (ATM_PDF->GetSKDetObj() == kDetObj_BeamSKDet) {
      std::cout << "-----------------------------" << std::endl;
      std::cout << "Setting Beam Detector covariance object.." << std::endl;
      ATM_PDF->setSkDetCov(SKDet);
    } else if (ATM_PDF->GetSKDetObj() == kDetObj_ATMPDDet) {
      std::cout << "-----------------------------" << std::endl;
      std::cout << "Setting ATMPD Detector covariance object.." << std::endl;
      ATM_PDF->setATMPDDetectorCov(ATMPDDet);
    }

    std::cout << "-----------------------------" << std::endl;
    std::cout << "Created Atm pdf with config:" << cfgName.Data() << std::endl;
    AtmPdfs.push_back(ATM_PDF);
    std::cout << "Atmospheric Sample: " << i << " | Name: " << ATM_PDF->GetSampleName() << " = OK" << std::endl;
    std::cout << "-------------------------------------------------------------------" << std::endl;
    nAtmSamplesLoaded += 1;

  }

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    if (AtmPdfs[iPDF]->GetSKDetObj() == kDetObj_AtmSKDet) {
      for (int PDFtoAdd=0;PDFtoAdd<AtmPdfs.size();PDFtoAdd++) {
	if (iPDF!=PDFtoAdd) {AtmPdfs[iPDF]->AddSamplePDF(AtmPdfs[PDFtoAdd]);}
      }
    }
  }

#if defined(CUDA)
    std::vector<samplePDFSKBaseGPU*> T2KPdfs;
#else
    std::vector<samplePDFSKBase*> T2KPdfs;
#endif

  if (useT2K) {
    std::cout << "-------------------------------------------------------------------" << std::endl;
    std::cout << "Loading T2K samples.." << "\n" << std::endl;

#if defined(CUDA)
    samplePDFSKBaseGPU *numu_pdf = new samplePDFSKBaseGPU(NuPOT, "configs/SKsamples_FHC1Rmu_summer2020.cfg", xsec);
    samplePDFSKBaseGPU *numubar_pdf = new samplePDFSKBaseGPU(NuBarPOT, "configs/SKsamples_RHC1Rmu_summer2020.cfg", xsec);
    samplePDFSKBaseGPU *nue_pdf = new samplePDFSKBaseGPU(NuPOT, "configs/SKsamples_FHC1Re_summer2020.cfg", xsec);
    samplePDFSKBaseGPU *nuebar_pdf = new samplePDFSKBaseGPU(NuBarPOT, "configs/SKsamples_RHC1Re_summer2020.cfg", xsec);
    samplePDFSKBaseGPU *nue1pi_pdf = new samplePDFSKBaseGPU(NuPOT, "configs/SKsamples_FHC1Re1de_summer2020.cfg", xsec);
    T2KPdfs.push_back(numu_pdf);
    T2KPdfs.push_back(numubar_pdf);
    T2KPdfs.push_back(nue_pdf);
    T2KPdfs.push_back(nuebar_pdf);
    T2KPdfs.push_back(nue1pi_pdf);
#else
    samplePDFSKBase *numu_pdf = new samplePDFSKBase(NuPOT, "configs/SKsamples_FHC1Rmu_summer2020.cfg", xsec);
    samplePDFSKBase *numubar_pdf = new samplePDFSKBase(NuBarPOT, "configs/SKsamples_RHC1Rmu_summer2020.cfg", xsec);
    samplePDFSKBase *nue_pdf = new samplePDFSKBase(NuPOT, "configs/SKsamples_FHC1Re_summer2020.cfg", xsec);
    samplePDFSKBase *nuebar_pdf = new samplePDFSKBase(NuBarPOT, "configs/SKsamples_RHC1Re_summer2020.cfg", xsec);
    samplePDFSKBase *nue1pi_pdf = new samplePDFSKBase(NuPOT, "configs/SKsamples_FHC1Re1de_summer2020.cfg", xsec);
    T2KPdfs.push_back(numu_pdf);
    T2KPdfs.push_back(numubar_pdf);
    T2KPdfs.push_back(nue_pdf);
    T2KPdfs.push_back(nuebar_pdf);
    T2KPdfs.push_back(nue1pi_pdf);
#endif

    for (int i=0;i<T2KPdfs.size();i++) {
      T2KPdfs[i]->setUseBeta(false);
      T2KPdfs[i]->setApplyBetaNue(false);
      T2KPdfs[i]->setApplyBetaDiag(false);
      T2KPdfs[i]->useNonDoubledAngles(true);
      std::cout << "-----------------------------" << std::endl;
      std::cout << "Setting Xsec covariance object.." << std::endl;
      T2KPdfs[i]->setSkXsecCov(xsec);

      if (T2KPdfs[i]->GetSKDetObj() == kDetObj_AtmSKDet) {
	std::cout << "-----------------------------" << std::endl;
	std::cout << "Setting AtmSKDet covariance object.." << std::endl;
	T2KPdfs[i]->setAtmSKDetCov(AtmDet);
	std::cout << "-----------------------------" << std::endl;
	std::cout << "Setting Supplementary T2K Detector covariance object.." << std::endl;
	T2KPdfs[i]->setSKDetBeamCov(SKDetBeam);
      } else if (T2KPdfs[i]->GetSKDetObj() == kDetObj_BeamSKDet) {
	std::cout << "-----------------------------" << std::endl;
	std::cout << "Setting Beam Detector covariance object.." << std::endl;
	T2KPdfs[i]->setSkDetCov(SKDet);
      } else if (T2KPdfs[i]->GetSKDetObj() == kDetObj_ATMPDDet) {
	std::cout << "-----------------------------" << std::endl;
	std::cout << "Setting ATMPD Detector covariance object.." << std::endl;
        T2KPdfs[i]->setATMPDDetectorCov(ATMPDDet);
      }

      std::cout << "-----------------------------" << std::endl;
      std::cout << "Beam Sample: " << i << " | Name: " << T2KPdfs[i]->GetSampleName() << " = OK" << std::endl;
      std::cout << "-------------------------------------------------------------------" << std::endl;
    }

    for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
      if (T2KPdfs[iPDF]->GetSKDetObj() == kDetObj_AtmSKDet) {
	for (int PDFtoAdd=0;PDFtoAdd<T2KPdfs.size();PDFtoAdd++) {
	  if (iPDF!=PDFtoAdd) {T2KPdfs[iPDF]->AddSamplePDF(T2KPdfs[PDFtoAdd]);}
	}
      }
    }

  }

  //####################################################################################
  std::cout << "-------------------------------------------------------------------" << std::endl;
  std::cout << "Moving onto analysis.. (For now this is just a test bed)" << "\n\n" << std::endl;

  std::vector<TString> SampleNames;
  std::vector<double> OscEventRate;
  std::vector<double> UnOscEventRate;
  std::vector<int> AtmSampleIDs(AtmPdfs.size());
  std::vector<double> AtmOscEventRates(AtmPdfs.size());

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    SampleNames.push_back(AtmPdfs[iPDF]->GetSampleName());
    AtmSampleIDs[iPDF] = AtmPdfs[iPDF]->GetSampleNumber();
  }
  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    SampleNames.push_back(T2KPdfs[iPDF]->GetSampleName());
  }

  std::cout << "-------------------------------------------------------------------" << std::endl;
  std::cout << "Reweighting to UnOscillated Spectra" << std::endl;

  osc->setParameters(oscpars_un);
  osc->acceptStep();
  osc->printNominalCurrProp();

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    AtmPdfs[iPDF]->reweight(osc->getPropPars());
  }

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    UnOscEventRate.push_back(AtmPdfs[iPDF]->get2DHist()->Integral());
  }

  PrintIntegral(AtmPdfs,WeightStyle,IncMigration,"./output/Atm_UnOscillatedEventRates.txt",Iden+"_UnOsc");

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    AtmPdfs[iPDF]->ResetSampleHistograms();
  }

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    T2KPdfs[iPDF]->reweight(osc->getPropPars());
  }

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    if (T2KPdfs[iPDF]->getNDim()==1) {
      UnOscEventRate.push_back(T2KPdfs[iPDF]->get1DHist()->Integral());
    } else {
      UnOscEventRate.push_back(T2KPdfs[iPDF]->get2DHist()->Integral());
    }
  }

  PrintIntegral(T2KPdfs,WeightStyle,IncMigration,"./output/T2K_UnOscillatedEventRates.txt",Iden+"_UnOsc");

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    T2KPdfs[iPDF]->ResetSampleHistograms();
  }

  std::cout << "-------------------------------------------------------------------" << std::endl;
  std::cout << "Reweighting to Oscillated Spectra" << std::endl;

  osc->setParameters(oscpars);
  osc->acceptStep();
  osc->printNominalCurrProp();

  for (int i=0;i<AtmPdfs.size();i++) {
    AtmPdfs[i]->reweight(osc->getPropPars());
  }

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    OscEventRate.push_back(AtmPdfs[iPDF]->get2DHist()->Integral());
  }

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    std::vector< std::vector<double> > Selection = AtmPdfs[iPDF]->GetStoredSelection();

    std::vector<double> SampleSelec(3);
    SampleSelec[0] = (double)kRecoSample;
    SampleSelec[1] = AtmPdfs[iPDF]->GetSampleNumber();
    SampleSelec[2] = AtmPdfs[iPDF]->GetSampleNumber()+1;
    Selection.push_back(SampleSelec);

    TH2D* Hist = AtmPdfs[iPDF]->get2DVarHist(kPDFBinning,kPDFBinning,Selection,false,WeightStyle);
    AtmOscEventRates[iPDF] = Hist->Integral();
  }
  SKCalibration->SetAtmOscEventRates(AtmSampleIDs, AtmOscEventRates);

  SKCalibration->PrintStoredEventRates();
  SKCalibration->CheckStoredEventRates();

  PrintIntegral(AtmPdfs,WeightStyle,IncMigration,"./output/Atm_OscillatedEventRates.txt",Iden+"_Osc");

  for (int iPDF=0;iPDF<AtmPdfs.size();iPDF++) {
    AtmPdfs[iPDF]->ResetSampleHistograms();
  }

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    T2KPdfs[iPDF]->reweight(osc->getPropPars());
  }

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    if (T2KPdfs[iPDF]->getNDim()==1) {
      OscEventRate.push_back(T2KPdfs[iPDF]->get1DHist()->Integral());
    } else {
      OscEventRate.push_back(T2KPdfs[iPDF]->get2DHist()->Integral());
    }
  }

  PrintIntegral(T2KPdfs,WeightStyle,IncMigration,"./output/T2K_OscillatedEventRates.txt",Iden+"_Osc");

  for (int iPDF=0;iPDF<T2KPdfs.size();iPDF++) {
    T2KPdfs[iPDF]->ResetSampleHistograms();
  }

  std::cout << "---------------------------------------------------------------" << std::endl;
  std::cout << std::setw(30) << "Sample Name" << " | " << std::setw(15) << "Oscillated" << " | " << std::setw(15) << "UnOscillated" << std::endl;
  std::cout << "---------------------------------------------------------------" << std::endl;
  for (int iSample=0;iSample<SampleNames.size();iSample++) {
    std::cout << std::setw(30) << SampleNames[iSample] << " | " << std::setw(15) << OscEventRate[iSample] << " | " << std::setw(15) << UnOscEventRate[iSample] << std::endl;
  }

  return 0;

}

#if defined(CUDA)
void PrintIntegral(std::vector<samplePDFSKBaseGPU*> PDFS, int WeightStyle, bool IncMigration, TString OutputName, TString Iden) {
#else
  void PrintIntegral(std::vector<samplePDFSKBase*> PDFS, int WeightStyle, bool IncMigration, TString OutputName, TString Iden) {
#endif

  std::vector<int> sampleNumbers(PDFS.size());
  for (int i=0;i<PDFS.size();i++) {
    sampleNumbers[i] = PDFS[i]->GetSampleNumber();
  }

  std::vector< std::vector<double> > EventRate;

  std::vector<TString> ModeNames(kMaCh3_nModes);
  for (int iMode=0;iMode<kMaCh3_nModes;iMode++) {
    ModeNames[iMode] = MaCh3mode_ToString((MaCh3_Mode)iMode);
  }

  std::vector<TString> OscChannelNames;

  for (int iPDF=0;iPDF<PDFS.size();iPDF++) {
    int sampleNumber = sampleNumbers[iPDF];

    int nOscChannels = PDFS[iPDF]->getNMCSamples();
    OscChannelNames = PDFS[iPDF]->getMCSampleNames();

    EventRate = std::vector< std::vector<double> >();
    EventRate.resize(kMaCh3_nModes);
    for (int iMode=0;iMode<kMaCh3_nModes;iMode++) {
      EventRate[iMode].resize(nOscChannels);
    }

    for (int iMode=0;iMode<kMaCh3_nModes;iMode++) {

      for (int iOscChannel=0;iOscChannel<nOscChannels;iOscChannel++) {

	std::vector< std::vector<double> > Selection = PDFS[iPDF]->GetStoredSelection();

	std::vector<double> Selec1(3);
	Selec1[0] = (double)kRecoSample;
	Selec1[1] = PDFS[iPDF]->GetSampleNumber();
	Selec1[2] = PDFS[iPDF]->GetSampleNumber()+1;

	Selection.push_back(Selec1);

	std::vector<double> Selec2(3);
	Selec2[0] = (double)kM3Mode;
	Selec2[1] = iMode;
	Selec2[2] = iMode+1;

	Selection.push_back(Selec2);

	std::vector<double> Selec3(3);
	Selec3[0] = (double)kOscChannel;
	Selec3[1] = iOscChannel;
	Selec3[2] = iOscChannel+1;

	Selection.push_back(Selec3);

        TH2D* Hist = (TH2D*)(PDFS[iPDF]->get2DVarHist(kPDFBinning,kPDFBinning,Selection,false))->Clone();;
        TAxis* XAxis = (TAxis*)Hist->GetXaxis()->Clone();
        TAxis* YAxis = (TAxis*)Hist->GetYaxis()->Clone();

        if (IncMigration) {
          for (int jPDF=0;jPDF<PDFS.size();jPDF++) {
            if (iPDF==jPDF) continue;
	    if (PDFS[jPDF]->GetSKDetObj() == kDetObj_AtmSKDet) {
	      Hist->Add((TH2D*)PDFS[jPDF]->get2DVarHist(kPDFBinning,kPDFBinning,Selection,false,WeightStyle,XAxis,YAxis)->Clone());
	    }
          }
        }

        double Integral = Hist->Integral();
        EventRate[iMode][iOscChannel] = Integral;
      }
    }

    PrintToScreen(PDFS[iPDF]->GetSampleName(),EventRate,OscChannelNames,ModeNames,Iden,OutputName);
  }

}

void PrintToScreen(TString SampleName,std::vector< std::vector<double> > IntegralList, std::vector<TString> OscChannelNames, std::vector<TString> ModeNames,TString Iden,TString OutputFileName) {
  int space = 14;

  int nModes = IntegralList.size();
  int nOscChannels = IntegralList[0].size();

  bool printToFile;
  if (OutputFileName.CompareTo("/dev/null")) {printToFile = true;}
  else {printToFile = false;}

  ofstream outfile;
  if (printToFile) {
    outfile.open(OutputFileName.Data(), std::ios_base::app);
    outfile.precision(7);
  }

  double PDFIntegral = 0;

  std::vector<double> ChannelIntegral;
  ChannelIntegral.resize(nOscChannels);
  for (unsigned int i=0;i<ChannelIntegral.size();i++) {ChannelIntegral[i] = 0.;}

  std::cout << "-------------------------------------------------" << std::endl;

  if (printToFile) {
    outfile << "\\begin{table}[ht]" << std::endl;
    outfile << "\\begin{center}" << std::endl;
    outfile << "\\caption{Integral breakdown for sample: " << SampleName << "}" << std::endl;
    outfile << "\\label{" << Iden << "_" << SampleName << "_EventRate}" << std::endl;
    TString nColumns;
    for (int i=0;i<nOscChannels;i++) {nColumns+="|c";}
    nColumns += "|c|";
    outfile << "\\begin{tabular}{|l" << nColumns.Data() << "}" << std::endl;
    outfile << "\\hline" << std::endl;
  }

  std::cout << "Integral breakdown for sample: " << SampleName << std::endl;
  std::cout << std::endl;

  if (printToFile) {outfile << std::setw(space) << "Mode:";}
  std::cout << std::setw(space) << "Mode:" << " ";
  for (int i=0;i<nOscChannels;i++) {
    if (printToFile) {outfile << "&" << std::setw(space) << MaCh3Utils::SKSampleName_toLatexString(OscChannelNames[i]) << " ";}
    std::cout << std::setw(space) << OscChannelNames[i] << " ";
  }
  if (printToFile) {outfile << "&" << std::setw(space) << "Total:" << "\\\\ \\hline" << std::endl;}
  std::cout << std::setw(space) << "Total:" << std::endl;

  for (unsigned int i=0;i<nModes;i++) {
    double ModeIntegral = 0;
    if (printToFile) {outfile << std::setw(space) << ModeNames[i];}
    std::cout << std::setw(space) << ModeNames[i] << " ";

    for (unsigned int j=0;j<IntegralList[i].size();j++) {
      double Integral = IntegralList[i][j];

      if (Integral<1e-100) {Integral=0;}

      ModeIntegral += Integral;
      ChannelIntegral[j] += Integral;
      PDFIntegral += Integral;

      if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",Integral) << " ";}
      std::cout << std::setw(space) << Integral << " ";
    }
    if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ModeIntegral) <<  " \\\\ \\hline" << std::endl;}
    std::cout << std::setw(space) << ModeIntegral << std::endl;
  }

  if (printToFile) {outfile << std::setw(space) << "Total:";}
  std::cout << std::setw(space) << "Total:" << " ";
  for (unsigned int i=0;i<ChannelIntegral.size();i++) {
    if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ChannelIntegral[i]) << " ";}
    std::cout << std::setw(space) << ChannelIntegral[i] << " ";
  }
  if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",PDFIntegral) << " \\\\ \\hline" << std::endl;}
  std::cout << std::setw(space) << PDFIntegral << std::endl;

  std::cout << "-------------------------------------------------" << std::endl;

  if (printToFile) {
    outfile << "\\end{tabular}" << std::endl;
    outfile << "\\end{center}" << std::endl;
    outfile << "\\end{table}" << std::endl;
  }

  std::cout << std::endl;

  if (printToFile) {
    outfile << std::endl;
    outfile.close();
  }

}