NCExtrapolationNS.cxx

Go to the documentation of this file.

//$Id: NCExtrapolationNS.cxx,v 1.40 2008/01/11 15:53:31 bckhouse Exp $
//
//NCExtrapolationNS.cxx
//
//class to hold pdfs for doing nccc separation
//
//B. Rebel 2/2007

#include "NCUtils/Extrapolation/NCExtrapolationNS.h"
#include "NCUtils/Extrapolation/NCEnergyBin.h"
#include "NCUtils/Cuts/NCAnalysisCutsNC.h"
#include "MessageService/MsgService.h"
#include "AnalysisNtuples/ANtpAnalysisInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"

#include "TH2.h"
#include "TSystem.h"
#include "TDirectory.h"
#include "TMath.h"
#include <cassert>

using namespace NCType;

// These three values correspond to shower_shift loops//
const int     n1     = 21;
const int     n0     = 0;
const double  les    = 0.80;
const double  hes    = 1.20;
const int     n1d    = 20;
const double size_e=(hes-les)/n1d;
// These three values correspond to nc shift loops//
const int     nc1    = 21;
const int     nc0    = 0;
const double  lnc    = 0.20;
const double  hnc    = 1.80;
const int     ncd    = 20;
const double size_nc=(hnc-lnc)/ncd;
// These three values correspond to normalization shift loops//
const int     nr1    = 11;//21;
const int     nr0    = 10;//0;
const double  lnr    = 0.95;
const double  hnr    = 1.05;
const int     nrd    = 20;
const double size_nr=(hnr-lnr)/nrd;

// These are the sigmas//
const double signorm = NCType::kParameterSigmas[NCType::kNormalization];
const double siges   = NCType::kParameterSigmas[NCType::kShowerEnergy];
const double signc   = NCType::kParameterSigmas[NCType::kNCBackground];

//far and near detector masses
const double kmass_f    = 52.97*486*0.0254*7.874+52.97*484*0.011*1.19879;
const double kmass_n_cc = 3.141592654*1.*1.*(1.710210*7.874+0.746800*1.20317); //reference to cc analysis
const double kmass_n    = 2.4930144*(1.710210*7.874+0.746800*1.20317)*(4.7368-1.728)/4.; //for nc analysis
const double kscalef    = 1.00;

CVSID("$Id: NCExtrapolationNS.cxx,v 1.40 2008/01/11 15:53:31 bckhouse Exp $");

//......................................................................
NCExtrapolationNS::NCExtrapolationNS() :
  NCExtrapolation::NCExtrapolation()
{
}

//.......................................................................
NCExtrapolationNS::NCExtrapolationNS(NCAnalysisUtils *utils,
                                     std::map<Int_t,double>& dataNear,
                                     std::map<Int_t,double>& mcNear,
                                     std::map<Int_t,double>& dataFar,
                                     std::map<Int_t,double>& mcFar,
                                     bool useCC,
                                     bool useNC) :
  NCExtrapolation(utils, dataNear, mcNear, dataFar, mcFar, useCC, useNC)
{

  fExtrapolationType = NCType::kExtrapolationNS;

  //set up the chi^2 histograms
  fchisqall   = new TH2F("chisqall","",
                         NCType::kNumUMu3SqrBins,
                         NCType::kUMu3SqrStart,
                         NCType::kUMu3SqrEnd,
                         NCType::kNumDeltaMSqrBins,
                         NCType::kDeltaMSqrStart,
                         NCType::kDeltaMSqrEnd);
  fchisqallnc = new TH2F("chisqallnc","",
                         NCType::kNumUS3SqrBins,
                         NCType::kUS3SqrStart,
                         NCType::kUS3SqrEnd,
                         NCType::kNumDeltaMSqrBins,
                         NCType::kDeltaMSqrStart,
                         NCType::kDeltaMSqrEnd);

}

//----------------------------------------------------------------------
NCExtrapolationNS::~NCExtrapolationNS()
{
}

//----------------------------------------------------------------------
//this method allows the user to pass in appropriate settings for the
//extrapolation such as locations of files, whether to generate the files
//again, etc.  the registry comes from the job module GetConfig method
void NCExtrapolationNS::Prepare(const Registry& r)
{
  MSG("NCExtrapolationNS", Msg::kWarning) << "NCExtrapolationNS::Prepare(const Registry& r)" << endl;
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string

  NCExtrapolation::Prepare(r);

  if(r.Get("Decay",                    tmpi)) fDecay                    = tmpi;
  if(r.Get("DoMC",                     tmpi)) fDoMC                     = tmpi;
  if(r.Get("BeamFileName",             tmps)) fBeamFileName             = tmps;
  if(r.Get("BeamFileNamePostShutDown", tmps)) fBeamFileNamePostShutDown = tmps;
  if(r.Get("TrueCorrectionHistNear",   tmps)) fTrueCorrectionHistNear   = tmps;
  if(r.Get("TrueCorrectionHistFar",    tmps)) fTrueCorrectionHistFar    = tmps;
  if(r.Get("FarTrueCorrectionFile",    tmps)) fFarTrueCorrectionFile    = tmps;
  if(r.Get("NearTrueCorrectionFile",   tmps)) fNearTrueCorrectionFile   = tmps;

  TString topDir="NCUtils/data";
  TString base="";
  base=getenv("SRT_PRIVATE_CONTEXT");
  if(base!="" && base!="."){
    // check if directory exists in SRT_PRIVATE_CONTEXT
    TString path = base + "/" + topDir;
    void *dir_ptr = gSystem->OpenDirectory(path);
    if(!dir_ptr) base=getenv("SRT_PUBLIC_CONTEXT");
  }
  // if it doesn't exist then use SRT_PUBLIC_CONTEXT
  else base=getenv("SRT_PUBLIC_CONTEXT");

  if(base==""){
    MSG("NCExtrapolationNS",Msg::kFatal) << "NCExtrapolationNS: No SRT_PUBLIC_CONTEXT set"
                                         << endl;
    assert(false);
  }
  topDir += "/";
  topDir.Prepend(base+"/");
  MSG("NCExtrapoationNS", Msg::kInfo)
    << "NCExtrapoationNS reading files from: " << topDir << endl;

  fBeamFileName.Prepend(topDir);
  fBeamFileNamePostShutDown.Prepend(topDir);
  fFarTrueCorrectionFile.Prepend(topDir);
  fNearTrueCorrectionFile.Prepend(topDir);

  MakeHistograms();

  tmpb = 0;
  tmpi = 0;
  tmpd = 0.;
  tmps = "";

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::MakeHistograms()
{
  fbinvect[0]=0.;
  fbinvect[1]=0.5;
  for(int i=2;i<200;i++){
    fbinvect[i]=0.25+(i)*0.25;
    if(i==199) fbinvect[i]=200.;
  }

  //save the current directory before opening some files
  TDirectory *savedir = gDirectory;

  TFile beamfile(fBeamFileName,"READ");
  fbmat_norm_r2 = dynamic_cast<TH2D *>(beamfile.Get("bmat_norm_r2"));
  fbmat_norm_r2->SetDirectory(savedir);
  TFile beamfile_post(fBeamFileNamePostShutDown,"READ");
  fbmat_norm_r2_post = dynamic_cast<TH2D *>(beamfile.Get("bmat_norm_r2"));
  fbmat_norm_r2_post->SetDirectory(savedir);

  savedir->cd();

  //these histograms will be made by the method now
  fSingleMCNearNuCC = new TH1D("singleMCNearNuCC", "", fnbins, fbinvect);
  fSingleMCFarNuCC = new TH1D("singleMCFarNuCC", "", fnbins, fbinvect);
  fSingleMCNearAllCC = new TH1D("singleMCNearAllCC", "", fnbins, fbinvect);
  fSingleMCFarAllCC = new TH1D("singleMCFarAllCC", "", fnbins, fbinvect);
  fMCNearNuCC = new TH1D("mcNearNuCC", "", fnbins, fbinvect);
  fMCFarNuCC = new TH1D("mcFarNuCC", "", fnbins, fbinvect);
  fMCNearAllCC = new TH1D("mcNearAllCC", "", fnbins, fbinvect);
  fMCFarAllCC = new TH1D("mcFarAllCC", "", fnbins, fbinvect);
  //LLH deprecating fcorf = new TH1D(fTrueCorrectionHistFar, "", fnbins, fbinvect); //for cc only - LLH
  //LLH deprecating fcor = new TH1D(fTrueCorrectionHistNear, "", fnbins, fbinvect); //for cc only - LLH
  fcorall = new TH1D("fcorall","",fnbins,fbinvect);  //for cc and nc? - LLH
  fcorfall = new TH1D("fcorfall","",fnbins,fbinvect);  //for cc and nc? - LLH

//   TFile corrf(fFarTrueCorrectionFile,"READ");
//   fcorf = dynamic_cast<TH1D *>(corrf.Get(fTrueCorrectionHistFar));
//   fcorf->SetDirectory(savedir);
//   TFile corr(fNearTrueCorrectionFile,"READ");
//   fcor = dynamic_cast<TH1D *>(corr.Get(fTrueCorrectionHistNear));
//   fcor->SetDirectory(savedir);

  fTrueND = new TH1F("trueND", "", fnbins, fbinvect);
  fTrueNDAfterRT = new TH1F("trueNDAfterRT", "", fnbins, fbinvect);
  fTrueNDAfterRT_noeffcorr = new TH1F("trueNDAfterRT_noeffcorr", "", fnbins, fbinvect); //LLH adding
  fRecoND = new TH1F("recoND", "", fnbins, fbinvect);
  fTrueFD = new TH1F("trueFD", "", fnbins, fbinvect);
  fTrueFD_nc = new TH1F("trueFD_nc", "", fnbins, fbinvect);
  fRecoFD = new TH1F("recoFD", "", fnbins, fbinvect);
  fRecoFD_nc = new TH1F("recoFD_nc", "", fnbins, fbinvect);
  fTrueFDAfterBM = new TH1F("trueFDAfterBM", "", fnbins, fbinvect);
  fTrueFDAfterBM_noacccorr = new TH1F("trueFDAfterBM_noacccorr", "", fnbins, fbinvect); //LLH adding
  fTrueFDAfterBM_posteffcorr = new TH1F("trueFDAfterBM_posteffcorr", "", fnbins, fbinvect); //LLH adding
  fTrueFDAfterOsc = new TH1F("trueFDAfterOsc", "", fnbins, fbinvect); //LLH adding
  fRecoFDAfterBM = new TH1F("recoFDAfterBM", "", fnbins, fbinvect);
//  fdatsel = new TH1F("datselstore", "", fnbins, fbinvect); //LLH adding, store before POT scaling

  //used in efficiency calculation in ND
  fallccmc = new TH1F("allccmc","",fnbins,fbinvect);
  fallccmcf = new TH1F("allccmcf","",fnbins,fbinvect); //LLH adding for proper eff calc
  fallncmcf = new TH1F("allncmcf","",fnbins,fbinvect); //LLH adding for proper eff calc

  // FAR CC AND NC TRUE-RECO
//not used-LLH ftruerecomat = new TH2F("truerecomat", "", fnbins, fbinvect, fnbins, fbinvect);
//not used-LLH ftruerecomatnc = new TH2F("truerecomatnc", "", fnbins, fbinvect, fnbins, fbinvect);
//not used-LLH ftruerecomatb = new TH2F("truerecomatb",  "", fnbins, fbinvect, fnbins, fbinvect);
//not used-LLH ftruerecomatncb = new TH2F("truerecomatncb","", fnbins, fbinvect, fnbins, fbinvect);

  // FAR CC AND NC SELECTION EFFICIENCIES (SIGNAL AND BACKGROUND)
  feffic = new TH1F("effic",   "",fnbins, fbinvect);
  fefficnc = new TH1F("efficnc", "",fnbins, fbinvect);
  fefficbg = new TH1F("efficb",  "",fnbins, fbinvect);
  fefficncbg = new TH1F("efficncb","",fnbins, fbinvect);
  fratccnc = new TH1F("ratccnc", "",fnbins, fbinvect);
  feffic_numcheck = new TH1F("effic_numcheck",   "",fnbins, fbinvect); //LLH adding to check eff calc

  // TAUS
  ftruerecomatt = new TH2F("truerecomatt", "", fnbins, fbinvect, fnbins, fbinvect);
  ftruerecomatbt = new TH2F("truerecomatbt", "", fnbins, fbinvect, fnbins, fbinvect);
  feffict = new TH1F("effict", "", fnbins, fbinvect);
  fefficbgt = new TH1F("efficbt", "", fnbins, fbinvect);
  frattau = new TH1F("rattau", "", fnbins, fbinvect);

  // NEAR
  frecotruemat            = new TH2F("recotruemat","",fnbins,fbinvect,fnbins,fbinvect);
  fpurityn                = new TH1F("purityn","",fnbins,fbinvect);
  fefficiencyn            = new TH1F("efficiencyn","",fnbins,fbinvect);

  // FAR CC AND NC TRUE-RECO  POST
//not used-LLH  ftruerecomat_post            = new TH2F("truerecomat_post",   "",fnbins,fbinvect,fnbins,fbinvect);
//not used-LLH  ftruerecomatnc_post          = new TH2F("truerecomatnc_post", "",fnbins,fbinvect,fnbins,fbinvect);
//not used-LLH  ftruerecomatb_post           = new TH2F("truerecomatb_post",  "",fnbins,fbinvect,fnbins,fbinvect);
//not used-LLH  ftruerecomatncb_post         = new TH2F("truerecomatncb_post","",fnbins,fbinvect,fnbins,fbinvect);

  // FAR CC AND NC SELECTION EFFICIENCIES (SIGNAL AND BACKGROUND)
  feffic_post                  = new TH1F("effic_post",   "",fnbins,fbinvect);
  fefficnc_post                = new TH1F("efficnc_post", "",fnbins,fbinvect);
  fefficbg_post                 = new TH1F("efficb_post",  "",fnbins,fbinvect);
  fefficncbg_post               = new TH1F("efficncb_post","",fnbins,fbinvect);
  fratccnc_post                = new TH1F("ratccnc_post", "",fnbins,fbinvect);
  // TAUS
  ftruerecomatt_post           = new TH2F("truerecomatt_post",   "",fnbins,fbinvect,fnbins,fbinvect);
  ftruerecomatbt_post          = new TH2F("truerecomatbt_post",  "",fnbins,fbinvect,fnbins,fbinvect);
  feffict_post                 = new TH1F("effict_post",   "",fnbins,fbinvect);
  fefficbgt_post                = new TH1F("efficbt_post",  "",fnbins,fbinvect);
  frattau_post                 = new TH1F("rattau_post","",fnbins,fbinvect);

  // Near
  frecotruemat_post            = new TH2F("recotruemat_post","",fnbins,fbinvect,fnbins,fbinvect);
  fpurityn_post                = new TH1F("purityn_post","",fnbins,fbinvect);
  fefficiencyn_post            = new TH1F("efficiencyn_post","",fnbins,fbinvect);

  //PREDICTED SPECTRUM CC - bestfit sys unosc
  //fnoosctruereal           = new TH1F("noosctruereal","",fnbins,fbinvect); //not being used - LLH
  fnoosctrue               = new TH1F("noosctrue","",fnbins,fbinvect); //unosc prediction cc, trueE
  fnooscreco               = new TH1F("nooscreco","",fnbins,fbinvect); //unosc prediction cc, recoE
  fnooscreconc             = new TH1F("nooscreconc","",fnbins,fbinvect); //unosc prediction nc bg, recoE
  fnooscrecotau            = new TH1F("nooscrecotau","",fnbins,fbinvect);

  //PREDICTED SPECTRUM NC - bestfit sys unosc
  fnoosctrue_ncfit             = new TH1F("noosctrue_ncfit","",fnbins,fbinvect); //unosc prediction nc, trueE, no sel eff!!!
  fnooscreco_ncfit            = new TH1F("nooscreco_ncfit","",fnbins,fbinvect); //unosc preditcion nc, recoE
  fnooscrecocc_ncfit             = new TH1F("nooscrecocc_ncfit","",fnbins,fbinvect); //onosc precition cc bg, recoE
  fnooscrecotau_ncfit          = new TH1F("nooscrecotau_ncfit","",fnbins,fbinvect);

  //PREDICTED SPECTRUM CC POST
  //fnoosctruereal_post          = new TH1F("noosctruereal_post","",fnbins,fbinvect);  //not being used - LLH
  fnoosctrue_post              = new TH1F("noosctrue_post","",fnbins,fbinvect);
  fnooscreco_post              = new TH1F("nooscreco_post","",fnbins,fbinvect);
  fnooscreconc_post            = new TH1F("nooscreconc_post","",fnbins,fbinvect);
  fnooscrecotau_post           = new TH1F("nooscrecotau_post","",fnbins,fbinvect);

  //PREDICTED SPECTRUM NC post
  fnoosctrue_ncfit_post            = new TH1F("noosctrue_ncfit_post","",fnbins,fbinvect);
  fnooscreco_ncfit_post           = new TH1F("nooscreco_ncfit_post","",fnbins,fbinvect);
  fnooscrecocc_ncfit_post            = new TH1F("nooscrecocc_ncfit_post","",fnbins,fbinvect);
  fnooscrecotau_ncfit_post         = new TH1F("nooscrecotau_ncfit_post","",fnbins,fbinvect);

  // auxilary spectra
  fnooscar              = new TH1F("nooscar","",fnbins,fbinvect);
  fnooscar_ccselcc      = new TH1F("nooscar_ccselcc","",fnbins,fbinvect);
  fnooscar_ccselnc         = new TH1F("nooscar_ccselnc","",fnbins,fbinvect);
  fnooscar_ncselnc      = new TH1F("nooscar_ncselnc","",fnbins,fbinvect);
  fnooscar_ncselcc         = new TH1F("nooscar_ncselcc","",fnbins,fbinvect);
  //Taus
  fnooscar_ccselcct        = new TH1F("nooscar_ccselcct","",fnbins,fbinvect);
  fnooscar_ccselnct     = new TH1F("nooscar_ccselnct","",fnbins,fbinvect);

  //POST
  fnooscar_post                = new TH1F("nooscar_post","",fnbins,fbinvect);
  fnooscar_ccselcc_post         = new TH1F("nooscar_ccselcc_post","",fnbins,fbinvect);
  fnooscar_ccselnc_post         = new TH1F("nooscar_ccselnc_post","",fnbins,fbinvect);
  fnooscar_ncselnc_post         = new TH1F("nooscar_ncselnc_post","",fnbins,fbinvect);
  fnooscar_ncselcc_post         = new TH1F("nooscar_ncselcc_post","",fnbins,fbinvect);
  //Taus
  fnooscar_ccselcct_post       = new TH1F("nooscar_ccselcct_post","",fnbins,fbinvect);
  fnooscar_ccselnct_post       = new TH1F("nooscar_ccselnct_post","",fnbins,fbinvect);

  // BEST FIT TRUE AND RECO FROM FIT  CC

  fbestfitreco            = new TH1F("bestfitreco","",fnbins,fbinvect);
  fbestfitreconc          = new TH1F("bestfitreconc","",fnbins,fbinvect);
  fbestfitrecotau         = new TH1F("bestfitrecotau","",fnbins,fbinvect);

  // BEST FIT TRUE AND RECO FROM FIT NC

  fbestfitreco_ncfit         = new TH1F("bestfitreco_ncfit","",fnbins,fbinvect);
  fbestfitrecocc_ncfit          = new TH1F("bestfitrecocc_ncfit","",fnbins,fbinvect);
  fbestfitrecotau_ncfit      = new TH1F("bestfitrecotau_ncfit","",fnbins,fbinvect);

  // BEST FIT TRUE AND RECO FROM FIT CC

  fbestfitreco_post         = new TH1F("bestfitreco_post","",fnbins,fbinvect);
  fbestfitreconc_post       = new TH1F("bestfitreconc_post","",fnbins,fbinvect);
  fbestfitrecotau_post      = new TH1F("bestfitrecotau_post","",fnbins,fbinvect);

  // BEST FIT TRUE AND RECO FROM FIT NC

  fbestfitreco_ncfit_post      = new TH1F("bestfitreco_ncfit_post","",fnbins,fbinvect);
  fbestfitrecocc_ncfit_post       = new TH1F("bestfitrecocc_ncfit_post","",fnbins,fbinvect);
  fbestfitrecotau_ncfit_post   = new TH1F("bestfitrecotau_ncfit_post","",fnbins,fbinvect);

  // HELP ON FITTING CC
  ftemptrue               = new TH1F("temptrue","",fnbins,fbinvect);
  ftempreco               = new TH1F("tempreco","",fnbins,fbinvect);
  ftempreconc             = new TH1F("tempreconc","",fnbins,fbinvect);
  ftemprecotau            = new TH1F("temprecotau","",fnbins,fbinvect);

  // HELP ON FITTING NC
  ftemptrue_ncfit            = new TH1F("temptrue_ncfit","",fnbins,fbinvect);
  ftempreco_ncfit            = new TH1F("tempreco_ncfit","",fnbins,fbinvect);
  ftemprecocc_ncfit             = new TH1F("temprecocc_ncfit","",fnbins,fbinvect);
  ftemprecotau_ncfit         = new TH1F("temprecotau_ncfit","",fnbins,fbinvect);

  // HELP ON FITTING CC
  ftemptrue_post            = new TH1F("temptrue_post","",fnbins,fbinvect);
  ftempreco_post            = new TH1F("tempreco_post","",fnbins,fbinvect);
  ftempreconc_post          = new TH1F("tempreconc_post","",fnbins,fbinvect);
  ftemprecotau_post         = new TH1F("temprecotau_post","",fnbins,fbinvect);

  // HELP ON FITTING NC
  ftemptrue_ncfit_post         = new TH1F("temptrue_ncfit_post","",fnbins,fbinvect);
  ftempreco_ncfit_post         = new TH1F("tempreco_ncfit_post","",fnbins,fbinvect);
  ftemprecocc_ncfit_post          = new TH1F("temprecocc_ncfit_post","",fnbins,fbinvect);
  ftemprecotau_ncfit_post      = new TH1F("temprecotau_ncfit_post","",fnbins,fbinvect);

  int nbins1 = 15;
  if(fDoMC) nbins1 = 39;

  float binvect1[40];
  binvect1[0]=0;
  binvect1[1]=1.;
  for (int i=2; i < nbins1+1; i++){
    if(fDoMC > 0 && i <= 19)  binvect1[i]=0.5+(i)*0.5;
    if(fDoMC > 0 && i > 19)   binvect1[i]=10.+(i-19)*1.;
    if(fDoMC < 1 && i <= 10)  binvect1[i]=(i)*1.;
    if(fDoMC < 1 && i > 10)   binvect1[i]=10.+(i-10)*2.;
  }

  fbestfitrecotau_rebinned        = new TH1F("bestfitrecotau_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreco_rebinned           = new TH1F("bestfitreco_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreconc_rebinned         = new TH1F("bestfitreconc_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  foscreco_rebinned                  = new TH1F("oscreco_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreco_rebinned              = new TH1F("tempreco_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreconc_rebinned            = new TH1F("tempreconc_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreco_rebinned             = new TH1F("nooscreco_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreconc_rebinned           = new TH1F("nooscreconc_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecotau_rebinned          = new TH1F("nooscrecotau_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);

  fbestfitrecotau_rebinned_post       = new TH1F("bestfitrecotau_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreco_rebinned_post          = new TH1F("bestfitreco_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreconc_rebinned_post        = new TH1F("bestfitreconc_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  foscreco_rebinned_post                     = new TH1F("oscreco_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreco_rebinned_post             = new TH1F("tempreco_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreconc_rebinned_post           = new TH1F("tempreconc_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreco_rebinned_post            = new TH1F("nooscreco_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreconc_rebinned_post          = new TH1F("nooscreconc_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecotau_rebinned_post         = new TH1F("nooscrecotau_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);

  //NC
  fbestfitrecotau_ncfit_rebinned     = new TH1F("bestfitrecotau_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreco_ncfit_rebinned        = new TH1F("bestfitreco_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitrecocc_ncfit_rebinned         = new TH1F("bestfitrecocc_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  foscreco_ncfit_rebinned            = new TH1F("oscreco_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreco_ncfit_rebinned           = new TH1F("tempreco_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftemprecocc_ncfit_rebinned            = new TH1F("temprecocc_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreco_ncfit_rebinned          = new TH1F("nooscreco_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecocc_ncfit_rebinned           = new TH1F("nooscrecocc_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecotau_ncfit_rebinned        = new TH1F("nooscrecotau_ncfit_rebinned","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);


  fbestfitrecotau_ncfit_rebinned_post    = new TH1F("bestfitrecotau_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitreco_ncfit_rebinned_post       = new TH1F("bestfitreco_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fbestfitrecocc_ncfit_rebinned_post        = new TH1F("bestfitrecovc_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  foscreco_ncfit_rebinned_post       = new TH1F("oscreco_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftempreco_ncfit_rebinned_post          = new TH1F("tempreco_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  ftemprecocc_ncfit_rebinned_post           = new TH1F("temprecocc_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscreco_ncfit_rebinned_post         = new TH1F("nooscreco_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecocc_ncfit_rebinned_post          = new TH1F("nooscrecocc_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
  fnooscrecotau_ncfit_rebinned_post       = new TH1F("nooscrecotau_ncfit_rebinned_post","",
                                        NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);

  return;
}

//----------------------------------------------------------------------
//add events to the histograms used for the weighting
void NCExtrapolationNS::AddEvent(ANtpHeaderInfo    *headerInfo,
                                 ANtpBeamInfo      *beamInfo,
                                 ANtpRecoInfo      *recoInfo,
                                 ANtpAnalysisInfo  *analysisInfo,
                                 ANtpTruthInfoBeam *truthInfo,
                                 int runType,
                                 bool useMCAsData,
                                 int fileType)
{

  //++++Energies are set to the correct NC/CC values in the NCExtrapolationModule
  //before being passed into the extrapolations

  //only use the L010z185i beam
  if(beamInfo->beamType != BeamType::ToZarko(BeamType::kL010z185i))
    return;

  if(headerInfo->detector == (int)Detector::kNear)
    AddNearEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                 truthInfo, runType, useMCAsData, fileType);

  if(headerInfo->detector == (int)Detector::kFar)
    AddFarEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                truthInfo, runType, useMCAsData, fileType);

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::AddNearEvent(ANtpHeaderInfo *headerInfo,
                                     ANtpBeamInfo *beamInfo,
                                     ANtpRecoInfo *recoInfo,
                                     ANtpAnalysisInfo *analysisInfo,
                                     ANtpTruthInfoBeam *truthInfo,
                                     int runType,
                                     bool useMCAsData,
                                     int fileType)
{

  //fill histogram for all MC events with a true vertex in the fiducial volume

  //Note: if this code is run on stripped uDST files, it will miss events that were reco'd
  //reco'd outside the FV!  Make sure to run on unstripped files or find a way to apply
  //a correction... This should be a small effect - LLH
  if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
    NCAnalysisCutsNC cuts;
    cuts.SetInfoObjects(headerInfo, beamInfo, 0, 0, 0, truthInfo);
    if( cuts.InFiducialVolumeTrue() && headerInfo->events != 0
        && truthInfo->interactionType == 1
        && TMath::Abs(truthInfo->nuFlavor) == 14){

      if(truthInfo->nuFlavor > 0) fallccmc->Fill(truthInfo->nuEnergy, recoInfo->weight);

      //check if this is a new snarl
//       if(headerInfo->newSnarl > 0) AddSingleMCEventsToHistogram(Detector::kNear);

      //Following lines were for fcorf calculation, since we can't do
      //this from uDST files, comment out for now - LLH

//       fTrueNuEnergy.push_back(truthInfo->nuEnergy);
//       fTrueNuFlavor.push_back(truthInfo->nuFlavor);
//       fTrueNuComplete.push_back(truthInfo->eventCompleteness);
    }//end if true vertex is in fiducial volume

  }//end if mc

  //Here apply FV and cleaning cuts!
  NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                            truthInfo, runType, useMCAsData, fileType);

  //dont bother with events outside the fiducial volume
//   if(recoInfo->inFiducialVolume < 1) return;

//   if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){

//     feshnd.push_back(recoInfo->showerEnergyCC);
//     femund.push_back(recoInfo->muEnergy);
//     fenutnd.push_back(truthInfo->nuEnergy);
//     fccnd.push_back(truthInfo->interactionType);
//     if(analysisInfo->isNC > 0)
//       fpidpnd.push_back(NCType::kNC);
//     else if(analysisInfo->isCC > 0)
//       fpidpnd.push_back(NCType::kCC);
//     fflavnd.push_back(truthInfo->nuFlavor);
//     finisnd.push_back(truthInfo->initialState);
//     fweind.push_back(recoInfo->weight);

//   }//end if MC

  return;
}

//----------------------------------------------------------------------
//all events passed to this method are garaunteed to be cc events, have
//their true interaction vertex within the fiducial volume and are either
//numu or numubar

//Note (9/10/07) This is deprecated for the moment since uDST's don't save
//info on generated MC events that weren't found by reconstruction - LLH
void NCExtrapolationNS::AddSingleMCEventsToHistogram(Detector::Detector_t det)
{
  //double loop over mc events in snarl and see if there are any repeats
  //for each entry with duplicates, store the number of duplicates found with j>i
  //store only instances where number of additional duplicates is 0
  //Because we are saving true energies, there's no need to worry about event completeness

  //note that if cleaning cut are applied to all events in the strip file,
  //then this method will not find any duplicates! - LLH
  vector<int> repeated;


  for(uint i = 0; i < fTrueNuEnergy.size(); ++i){
    repeated.push_back(-1);
    uint duplicatesCtr = 0; //for each duplicate found for entry i, increment this by 1

    for(uint j = i+1; j < fTrueNuEnergy.size(); ++j) //only search through j > i
    {
       if(fTrueNuEnergy[i] == fTrueNuEnergy[j])
          duplicatesCtr++;

    }//end inner loop

    repeated[i] = duplicatesCtr; //storing number of duplicates found after entry i

    MSG("NCExtrapolationNS", Msg::kInfo) << "i= " << i  //LLH changed from kDebug, 100
                                         << ", repeated[i] = " << repeated[i]
                                         << ", duplicatesCtr = " << duplicatesCtr
                                         <<", energy = " << fTrueNuEnergy[i]
                                         << endl;

  }//end outer loop

  //loop over the array one more time and fill the histogram
  //with single events or the event with the largest completeness
  //if it is a duplicate
  for(uint i = 0; i < fTrueNuEnergy.size(); ++i){

    if(det == Detector::kNear){
       if(repeated[i] ==0) fSingleMCNearAllCC->Fill(fTrueNuEnergy[i]);  //all cc-like
       fMCNearAllCC->Fill(fTrueNuEnergy[i]);
       if(fTrueNuFlavor[i] == 14){
          if(repeated[i] == 0) fSingleMCNearNuCC->Fill(fTrueNuEnergy[i]);  //true numu cc-like
          fMCNearNuCC->Fill(fTrueNuEnergy[i]);
      }//end if numu
    }//end if near

    else if(det == Detector::kFar){
      if(repeated[i] ==0) fSingleMCFarAllCC->Fill(fTrueNuEnergy[i]);  //all cc-like
      fSingleMCFarAllCC->Fill(fTrueNuEnergy[i]);
      if(fTrueNuFlavor[i] == 14){
        if(repeated[i] == 0) fSingleMCFarNuCC->Fill(fTrueNuEnergy[i]);  //true numu cc-like
        fMCFarNuCC->Fill(fTrueNuEnergy[i]);
      }//end if numu
    }//end if far
  }//end loop over events

  //reset the arrays (for next snarl)
  fTrueNuEnergy.clear();
  fTrueNuFlavor.clear();
  fTrueNuComplete.clear();

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::AddFarEvent(ANtpHeaderInfo *headerInfo,
                                    ANtpBeamInfo *beamInfo,
                                    ANtpRecoInfo *recoInfo,
                                    ANtpAnalysisInfo *analysisInfo,
                                    ANtpTruthInfoBeam *truthInfo,
                                    int runType,
                                    bool useMCAsData,
                                    int fileType)
{
   MAXMSG("NCExtrapolationNS", Msg::kInfo, 1) << "AddFarEvent " << endl;

  //fill histogram for all MC events with a true vertex in the fiducial volume
   if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){

      //don't double count and make sure event passed reco because
      //reconstruction efficiency is accounted for in the corf function
      //(eventually take this out)
      if(headerInfo->newSnarl && headerInfo->events!=0)
      {
         //Far detector efficiency corrections must be wrt all (numu) events in far detector
         if(truthInfo->nuFlavor == 14)
         {
            if(truthInfo->interactionType == 1) fallccmcf->Fill(truthInfo->nuEnergy, recoInfo->weight);

            //Following lines were for fcorf calculation, since we can't do
            //this from uDST files, comment out for now - LLH

            //check if this is a new snarl
//      if(headerInfo->newSnarl > 0) AddSingleMCEventsToHistogram(Detector::kFar);
//       fTrueNuEnergy.push_back(truthInfo->nuEnergy);
//       fTrueNuFlavor.push_back(truthInfo->nuFlavor);
//       fTrueNuComplete.push_back(truthInfo->eventCompleteness);
         }//end if true vertex is a numu

         //For NC events store all nu types
         if(truthInfo->interactionType == 0) fallncmcf->Fill(truthInfo->nuEnergy, recoInfo->weight);
      } //end if nonzero events and is new snarl

  }//end if mc


  NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                            truthInfo, runType, useMCAsData, fileType);

  return;
}

//----------------------------------------------------------------------
//just have to multiply total MC events with true vertex in the
//fiducial volume by the correction factor for split events
void NCExtrapolationNS::PrepareNearDetector()
{
  //get the correction histograms for the single events
//  fcor->Divide(fSingleMCNearNuCC, fMCNearNuCC);
//  fcorall->Divide(fSingleMCNearAllCC, fMCNearAllCC);

   cout <<"In PrepareNearDetector" << endl;

  //For now read in fcor and fcorall from Niki's files, since we don't have the
  //capability to do this from uDST's
//  TFile *corr      = new TFile("$MINOS_DATA/d84/niki/mad_cedar_daikon/near/correction_plots/near/norm_fix_near.root","READ");
   TFile *corr      = new TFile("/afs/fnal.gov/files/data/minos/d79/llhsu/eventFindingEff/norm_fix_near_uDST_strip.root","READ"); //using my own file - LLH with zero event cut
  fcor = (TH1F*)(corr->Get("cor_neu"));

  //Correct the true CC events for reconstruction efficiencies and event splitting
  fallccmc->Multiply(fcor);
  return;
}

//----------------------------------------------------------------------
//this is where you call the guts of the code to do the prediction and
//fit
void NCExtrapolationNS::PredictSpectrum()
{

   cout <<"In NCExtrapolationNS::PredictSpectrum!" << endl;

  //make the correction histograms used by the method
//LLH deprecating  fcorf->Divide(fSingleMCFarNuCC, fMCFarNuCC);
//LLH deprecating  fcorfall->Divide(fSingleMCFarAllCC, fMCFarAllCC);
  TFile *corrf      = new TFile("$MINOS_DATA/d84/niki/mad_cedar_daikon/near/correction_plots/far/norm_fix_far.root","READ");
  fcorf = (TH1F*)(corrf->Get("cor_neuf"));

  //find the position in the beam vector
  int runI = NCType::kRunI;
  int runII = runI;//NCType::kRunII; Temp! FIXME - LLH
  int farBeamI = NCExtrapolation::FindBeamIndex(Detector::kFar, BeamType::kL010z185i, runI);
  int farBeamII = NCExtrapolation::FindBeamIndex(Detector::kFar, BeamType::kL010z185i, runII);
  int nearBeamI = NCExtrapolation::FindBeamIndex(Detector::kNear, BeamType::kL010z185i, runI);
  int nearBeamII = NCExtrapolation::FindBeamIndex(Detector::kNear, BeamType::kL010z185i, runII);
  if(farBeamI > 100 || farBeamII > 100
     || nearBeamI > 100 || nearBeamII > 100){
    MSG("NCExtrapolationNS", Msg::kWarning) << "could not find "
                                            << "FD/ND beam in vector - bail"
                                            << " far: " << farBeamI << " " << farBeamII
                                            << " near: " << nearBeamI << " " << nearBeamII
                                            << endl;
    return;
  }

  //the next 2 lines just make everything look good without doing any fits
  //ie just the nominal mc and data
  //return;

  //get near and far histograms
  //datsel histograms are near data
  TH1F *datsel = new TH1F("datsel", "", fnbins, fbinvect);
  TH1F *datselpost = new TH1F("datselpost", "", fnbins, fbinvect);
  fdatselFD = new TH1F("datselstoreFD", "", fnbins, fbinvect); //LLH added
  fdatselFD_nc = new TH1F("datselstoreFD_nc", "", fnbins, fbinvect); //LLH added
  //TH1F *datselnc = fBeams[nearBeamI].GetDataHistogram(NCType::kNC);
  //TH1F *datselncpost = fBeams[nearBeamII].GetDataHistogram(NCType::kNC);
  FillDataHistogramFromEnergyBins(nearBeamI, NCType::kCC, datsel);
  FillDataHistogramFromEnergyBins(nearBeamII, NCType::kCC, datselpost);
  FillDataHistogramFromEnergyBins(farBeamI, NCType::kCC, fdatselFD); //LLH added
  FillDataHistogramFromEnergyBins(farBeamI, NCType::kNC, fdatselFD_nc); //LLH added

  //put far data histograms into Niki's names for the fitting
  FillDataHistogramFromEnergyBins(farBeamI, NCType::kCC, foscreco_rebinned);
  FillDataHistogramFromEnergyBins(farBeamII, NCType::kCC, foscreco_rebinned_post);
  FillDataHistogramFromEnergyBins(farBeamI, NCType::kNC, foscreco_ncfit_rebinned);
  FillDataHistogramFromEnergyBins(farBeamII, NCType::kNC, foscreco_ncfit_rebinned_post);

  //this method needs to know the relative scaling between near and far
  //exposures.  the matrix extrapolates based on the ND exposure
  //scale near detector exposure to be same as far detector
  double potNear = fBeams[nearBeamI].GetDataPOT();
  double potFar = fBeams[farBeamI].GetDataPOT();

  fNearExposureScale = potFar/potNear;
  cout <<"Scaling datsel!!!!!!! scale factor = " << fNearExposureScale
       <<"\npotFar = " << potFar
       <<"\npotNear = " << potNear
       << endl; //LLH
  fdatselND = (TH1F*)datsel->Clone("datselstoreND"); //LLH unmodified datsel
  datsel->Scale(fNearExposureScale);
  datselpost->Scale(fNearExposureScale);

  double scale_factNC = 1.;
  double nc_shift  = 1.;
  double enu_shift = 1.;

  // fit loop
  float dmsqmin     = 0.;
  float s2tmin      = 0.;
  float s2smin      = 0.;
  float chisqmin    = 9999.;
  float norm_min    = 9999.;
  float enu_shiftmin= 9999.;
  float nc_shiftmin = 9999.;
  int best_enu      = 9999;
  int best_nc       = 9999;

  int ndf = 39 + 39 - 3;
  if(fDoMC < 1) ndf=15+15-3;

  vector <vector <vector <double> > > chisqarraymins;
  chisqarraymins.resize(NCType::kNumDeltaMSqrBins);
  for(int i = 0; i < NCType::kNumDeltaMSqrBins; i++){
    chisqarraymins[i].resize(NCType::kNumUS3SqrBins);
    for(int j = 0; j < NCType::kNumUS3SqrBins; j++){
      chisqarraymins[i][j].resize(NCType::kNumUMu3SqrBins);
      for(int k = 0; k < NCType::kNumUMu3SqrBins; k++){
        chisqarraymins[i][j][k] = 1.e12;
      }
    }
  }

  float dmsqfit = 0.;
  float s2tfit = 0.;
  float s2sfit = 0.;

  TString math;

  TString noosc;

  noosctruearray.resize(n1);
  noosctruearraypost.resize(n1);

  //initializing arrays
  for(int i = 0; i < n1; ++i){
    math = TString::Format("truerecomat%d",i);
    matarray.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatnc%d",i);
    matarraync.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatb%d",i);
    matarrayb.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatncb%d",i);
    matarrayncb.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    //TAUS
    math = TString::Format("truerecomatt%d",i);
    matarrayt.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatbt%d",i);
    matarraybt.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));

    math = TString::Format("truerecomat_post%d",i);
    matarray_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatnc_post%d",i);
    matarraync_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatb_post%d",i);
    matarrayb_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatncb_post%d",i);
    matarrayncb_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    //TAUS
    math = TString::Format("truerecomatt_post%d",i);
    matarrayt_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
    math = TString::Format("truerecomatbt_post%d",i);
    matarraybt_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));

    for(int j = 0; j < nc1; ++j){

      noosc = TString::Format("noosctrue_%d_%d",i,j);
      noosctruearray[i].push_back(new TH1F(noosc,"",fnbins,fbinvect));
      noosc = TString::Format("noosctrue_post_%d_%d",i,j);
      noosctruearraypost[i].push_back(new TH1F(noosc,"",fnbins,fbinvect));
    }
  }//end initialization of arrays

  //Getting efficiencies (Far only)
  GetEff(feffic, fefficnc, fefficbg, fefficncbg, fratccnc, frattau, farBeamI);
  GetEff(feffic_post, fefficnc_post, fefficbg_post, fefficncbg_post, fratccnc_post, frattau_post, farBeamII);
  if(fDoMC==0){
    GetEfft(feffict, fefficbgt, farBeamI);
    GetEfft(feffict_post, fefficbgt_post, farBeamII);
  }

  double purbin = 0.;
  double newbin = 0.;

  // loop for truereco mat efficiency and purity For Far and Near
  int start = n1d/2;
  int end = start + 1;
  int startnc = ncd/2;
  int endnc = startnc + 1;

  if(fSystematicsToUse[NCType::kShowerEnergy]){
    start   = n0;
    end     = n1;
  }
  if(fSystematicsToUse[NCType::kNCBackground]){
    startnc = nc0;
    endnc   = nc1;
  }

  //use TH1::Copy function to copy the ftruerecomat* matricies into the
  //array of matricies matarray*

  for(int ii = start; ii < end; ii++){
    // SHOWER LOOP
    enu_shift = les+ii*size_e;
    MSG("NCExtrapolationNS", Msg::kInfo) << " enu  shift " << enu_shift << endl;
    cout <<"First Shower Loop!" << endl;

    RecoTrueEffPur(matarray[ii], matarraync[ii], matarrayb[ii],
                   matarrayncb[ii], enu_shift, farBeamI);
    // here TAUS
    if(fDoMC==0)
      RecoTrueEffPurTau(matarrayt[ii], matarraybt[ii], enu_shift, farBeamI);

    RecoTruePur(frecotruemat, fpurityn, fefficiencyn, enu_shift, nearBeamI);

    for(int gg = 0; gg < fnbins; gg++){
      purbin = fpurityn->GetBinContent(gg+1);
      newbin = purbin/(purbin+scale_factNC*(1-purbin));
      fpurityn->SetBinContent(gg+1,newbin);
    }

    // post shutdown
//     RecoTrueEffPur(matarray_post[ii], matarraync_post[ii], matarrayb_post[ii],
//                 matarrayncb_post[ii], enu_shift, farBeamII);

//     // here TAUS
//     if(fDoMC==0)
//       RecoTrueEffPurTau(matarrayt_post[ii], matarraybt_post[ii], enu_shift, farBeamII);

//     RecoTruePur(frecotruemat_post, f_purityn_post, fefficiencyn_post, enu_shift, nearBeamII);

//     for(int gg = 0; gg < fnbins; gg++){
//       purbin = fpurityn_post->GetBinContent(gg+1);
//       newbin = purbin/(purbin+scale_factNC*(1-purbin));
//       fpurityn_post->SetBinContent(gg+1,newbin);
//     }

    // NC LOOP
    for(int jj = startnc; jj < endnc; jj++){
      nc_shift = lnc+jj*size_nc;
      MSG("NCExtrapolationNS", Msg::kInfo) << " nc  shift " << nc_shift << endl;
      cout <<"NC loop within first shower loop" << endl;

      // MAKE PREDICTION UNOSCILLATED //
      for(int gg = 0; gg < fnbins; gg++){
        purbin = fpurityn->GetBinContent(gg+1);
        newbin = purbin/(purbin+nc_shift*(1-purbin));
        fpurityn->SetBinContent(gg+1,newbin);
      }

      //Temp - Code speedup for debugging only!
      Int_t bypass = 1;
      if(bypass == 1)
      {
         cout <<"Bypass 3!" << endl;
         TFile* f = new TFile("/local/scratch04/llhsu/tempFiles/L010z185i/extrapolation_fn_NS_16kND488FD_ncfvcor.root","READ");
         TString noosctruename = TString::Format("noosctrue_%d_%d", ii, jj);
         TH1F* tempHisto = (TH1F*)f->Get("plotsBeamMatrix/"+noosctruename);
         noosctruename += "_clone";
         noosctruearray[ii][jj] = (TH1F*)tempHisto->Clone(noosctruename);
         noosctruearray[ii][jj]->SetDirectory(0);
         //       fdatselND = (TH1F*)((f->Get("datselstoreND"))->Clone("datselstoreND"));
//         fpurityn = (TH1F*)((f->Get("purityn"))->Clone("purityn"));
//         fefficiencyn = (TH1F*)((f->Get("efficiencyn"))->Clone("efficiencyn"));
//         frecotruemat = (TH2F*)((f->Get("recotruemat"))->Clone("recotruemat"));
         //f->Close();
         cout <<"Done 3!" << endl;
      }
      else
      {

      //take ND "data" and generate predicted true far spectrum
      MakePrediction(datsel, fpurityn, fefficiencyn,
                     frecotruemat, noosctruearray[ii][jj], runI);

     } //remove else brackets with bypass

      // MAKE PREDICTION UNOSCILLATED //
//       for(int gg = 0; gg < fnbins; gg++){
//      purbin = fpuritynp->GetBinContent(gg+1);
//      newbin = purbin/(purbin+nc_shift*(1-purbin));
//      fpuritynp->SetBinContent(gg+1,newbin);
//       }

//       MakePrediction(datselpost, fpurityn_post, fefficiencyn_post,
//                   frecotruemat_post, noosctruearraypost[ii][jj], runII);

      // restore purity
      for(int gg = 0; gg < fnbins; gg++){
        purbin = fpurityn->GetBinContent(gg+1);
        newbin = (nc_shift*purbin)/(1-purbin+nc_shift*purbin);
        fpurityn->SetBinContent(gg+1,newbin);

//      purbin = fpuritynp->GetBinContent(gg+1);
//      newbin = (nc_shift*purbin)/(1-purbin+nc_shift*purbin);
//      fpuritynp->SetBinContent(gg+1,newbin);
      }
    }//end nc loop

//Temp LLH, for 1 iteration only!    frecotruemat->Reset();
//Temp LLH, for 1 iteration only!    fpurityn->Reset();
//Temp LLH, for 1 iteration only!    fefficiencyn->Reset();

//     frecotruemat_post->Reset();
//     fpurityn_post->Reset();
//     fefficiencyn_post->Reset();

  } // end shower loop

  MSG("NCExtrapolationNS", Msg::kInfo) << " DONE EFF PUR RECO->TRUE NEAR FAR " << endl;

  float chisqtemp1 = 1.e12;
  //  float chisqtemp2 = 1.e12;
  float chisqtemp3 = 1.e12;
  //  float chisqtemp4 = 1.e12;
  float chisqmin_norm = 1.e12;

  float nc_min   = 9999.;
  float sh_min   = 9999.;
  double norm = 1.;
  int decayMode = 1;
  int decayModeNC = 0;
  int decayModeTau = 11;
  if(fDecay == 1){
    decayMode = -1;
    decayModeNC = -1;
    decayModeTau = -11;
  }

  MSG("NCExtrapolationNS", Msg::kInfo) << " outside loop SHOWER Start "
                                       << start << " END " << end
                                       << " NC Start " << startnc
                                       << " END " << endnc << endl;


  int maxsin = 1;//99;//NCType::kNumUMu3SqrBins;
  int maxdm = 1; //119;//NCType::kNumDeltaMSqrBins;
  int maxfs = 1;//NCType::kNumUS3SqrBins;
  vector<double> sinSqrDeltaMSqr;

  //Fitting and oscillations??
  // SHOWER LOOP // ii
  for(int ii = start; ii < end; ii++){
    enu_shift = les + ii*size_e;
    MSG("NCExtrapolationNS", Msg::kDebug) << " Shower Loop " << ii << endl;
    cout <<"starting second shower loop, enu_shift = " << enu_shift << endl;
    // NC LOOP  // jj
    for(int jj = startnc; jj < endnc; jj++){
      nc_shift = lnc + jj*size_nc;
      MSG("NCExtrapolationNS", Msg::kDebug) << " NC Loop " << jj
                                            << " nc_shift " << nc_shift << endl;
      cout  <<"nc loop within second shower loop, nc_shift = " << nc_shift << endl;

      //I need four histos here
      //true cc sel as cc, true cc sel as nc,
      //true nc sel as nc, true nc sel as cc

      fnooscar->Reset();
      fnooscar_post->Reset();
      fnooscar_ccselcc->Reset();
      fnooscar_ccselcc_post->Reset();
      fnooscar_ccselnc->Reset();
      fnooscar_ccselnc_post->Reset();
      fnooscar_ncselnc->Reset();
      fnooscar_ncselnc_post->Reset();
      fnooscar_ncselcc->Reset();
      fnooscar_ncselcc_post->Reset();
      // taus
      fnooscar_ccselcct->Reset();
      fnooscar_ccselcct_post->Reset();
      fnooscar_ccselnct->Reset();
      fnooscar_ccselnct_post->Reset();

//Temp disabling to check fitting      fnooscar->Add(noosctruearray[ii][jj]); //Copying prediction into working histo!!
      fnooscar->Add(fallccmcf); //Temp using true FD spectrum (before pid efficiency corr) instead of prediction!

      fTrueFDAfterBM->Add(noosctruearray[ii][jj]);  //for storage
//       fnooscar_post->Add(noosctruearraypost[ii][jj]);

      //PRE
      fnooscar_ccselcc->Multiply(fnooscar,feffic);  //truecc, cclike, efficiency correction
      fTrueFDAfterBM_posteffcorr=(TH1F*)fnooscar_ccselcc->Clone("trueFDAfterBM_posteffcorr"); //LLH adding to compare to fTrueFD after PID'as applied
      fnooscar_ccselnc->Multiply(fnooscar,fefficbg); //truecc, nclike, efficiency correction
      //Obtain true NC Spectrum from the NC/CC Ratio
      fnooscar_ncselnc->Multiply(fnooscar,fratccnc); //truenc, nclike, cross-section correction
      fnooscar_ncselcc->Multiply(fnooscar,fratccnc); //truenc, cclike, cross-section correction
      fnooscar_ncselnc->Multiply(fefficnc); //efficiency correction
      fnooscar_ncselcc->Multiply(fefficncbg); //efficiency correction

      // POST
//       fnooscar_ccselcc_post->Multiply(fnooscar_post,fefficp);
//       fnooscar_ccselnc_post->Multiply(fnooscar_post,fefficbgp);
//       //Obtain true NC Spectrum from the NC/CC Ratio
//       fnooscar_ncselnc_post->Multiply(fnooscar_post,fratccncp);
//       fnooscar_ncselcc_post->Multiply(fnooscar_post,fratccncp);
//       fnooscar_ncselnc_post->Multiply(fefficncp);
//       fnooscar_ncselcc_post->Multiply(fefficncbgp);

      // taus
      if(fDoMC==0){
        // PRE
        fnooscar_ccselcct->Multiply(fnooscar,frattau);
        fnooscar_ccselcct->Multiply(feffict);
        fnooscar_ccselnct->Multiply(fnooscar,frattau);
        fnooscar_ccselnct->Multiply(fefficbgt);

        // POST
//      fnooscar_ccselcct_post->Multiply(fnooscar_post,frattaup);
//      fnooscar_ccselcct_post->Multiply(feffictp);
//      fnooscar_ccselnct_post->Multiply(fnooscar_post,frattaup);
//      fnooscar_ccselnct_post->Multiply(fefficbgtp);
      }

      //Loop over deltamsq
      for(int i = 0; i < maxdm; i++){
        dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;
        dmsqfit *= 1.e-3;

        MSG("NCExtrapolationNS", Msg::kInfo) << " on Delta m^2 = " << dmsqfit
                                             << " chi^2 min = " << chisqmin << endl;
        cout <<"Inside deltamsq loop! dmsqfit = " << dmsqfit << endl;

        sinSqrDeltaMSqr.clear();
        //really just need the value of the sin^2(1.27*\Delta m^2 L/E) factor
        //at the middle of each bin as the bins are of size 0.25 GeV in the
        //region of oscillation interest - that is good enough to get the oscillated
        //true value to within 0.0004%
        for(int dmb = 0; dmb < fnbins; ++dmb){  //Why was this set to 27?? - LLH
          sinSqrDeltaMSqr.push_back(TMath::Power(TMath::Sin(NCType::k127*dmsqfit
                                                            *NCType::kBaseLineFar/fnooscreco->GetBinCenter(dmb+1)),2.));
        }//end loop over energy bins

        for(int j = 0; j < maxsin; j++){
          MSG("NCExtrapolationNS", Msg::kDebug) << " In sin loop " << j <<  endl;
          s2tfit = (1.*j)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart; //Temp! (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;
          cout <<"Inside sinsq loop! s2tfit = " << s2tfit << endl;

          for(int k = 0; k < maxfs; k++){
            s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;
            cout <<"Inside fsterile loop! s2sfit = " << s2sfit << endl;

            //PRE SHUTDOWN

            //Oscillate true cc's selected as cc's after efficiency and turn them to reco
            // LAST 1 IS FOR CC
            //oscweight(fnooscar_ccselcc,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcc,ftemptrue,14,NCType::kCC);
            fTrueFDAfterOsc = (TH1F*)ftemptrue->Clone("trueFDAfterOsc"); //should be identical to fTrueFDAfterBM_posteffcorr w/o osc
            //truetoreco(fTrueFD,*ftempreco,matarray[ii]);  //Temp! To debug unsmearing - LLH
            truetoreco(ftemptrue,*ftempreco,matarray[ii]);
            fRecoFDAfterBM->Add(ftempreco);
            ftemptrue->Reset();
            //Oscillate true nc's selected as cc's after efficiency and turn them to reco
            //oscweight (fnooscar_ncselcc,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselcc,ftemptrue,14,NCType::kNC);
            truetoreco(ftemptrue,*ftempreconc,matarrayncb[ii]);
            ftempreconc->Scale(nc_shift);
            ftemptrue->Reset();
            //Oscillate true nc's selected as nc's after efficiency and turn them to reco
            //oscweight(fnooscar_ncselnc,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselnc,ftemptrue,14,NCType::kNC);
            truetoreco(ftemptrue,*ftempreco_ncfit,matarraync[ii]);
            ftempreco_ncfit->Scale(nc_shift);
            ftemptrue->Reset();
            //Oscillate true cc's selected as nc's after efficiency and turn them to reco
            //oscweight(fnooscar_ccselnc,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnc,ftemptrue,14,NCType::kCC);
            truetoreco(ftemptrue,*ftemprecocc_ncfit,matarrayb[ii]);
            ftemptrue->Reset();

            // TAUS  HERE
            if(fDoMC==0){
              //Oscillate true cc's selected as cc's after efficiency and turn them to reco
              // LAST 11 IS FOR CC TAU
              //oscweight(fnooscar_ccselcct,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
              oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcct,ftemptrue,16,NCType::kCC);
              truetoreco(ftemptrue,*ftemprecotau,matarrayt[ii]);
              ftemptrue->Reset();

              //Oscillate true cc's selected as nc's after efficiency and turn them to reco
              //oscweight(fnooscar_ccselnct,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
              oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnct,ftemptrue,16,NCType::kNC);
              truetoreco(ftemptrue,*ftemprecotau_ncfit,matarraybt[ii]);
              ftemptrue->Reset();
            }   //end taus

//          //POST SHUTDOWN

//          //Oscillate true cc's selected as cc's after efficiency and turn them to reco
//          // LAST 1 IS FOR CC
//          //oscweight(fnooscar_ccselcc_post,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
//          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcc_post,ftemptrue,14,NCType::kCC);
//          truetoreco(ftemptrue,*ftempreco_post,matarray_post[ii]);
//          ftemptrue->Reset();
//          //Oscillate true nc's selected as cc's after efficiency and turn them to reco
//          //oscweight(fnooscar_ncselcc_post,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
//          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselcc_post,ftemptrue,14,NCType::kNC);
//          truetoreco(ftemptrue,*ftempreconc_post,matarrayncb_post[ii]);
//          ftempreconc_post->Scale(nc_shift);
//          ftemptrue->Reset();

//          //Oscillate true nc's selected as nc's after efficiency and turn them to reco
//          //oscweight(fnooscar_ncselnc_post,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
//          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselnc_post,ftemptrue,14,NCType::kNC);
//          truetoreco(ftemptrue,*ftempreco_ncfit_post,matarraync_post[ii]);
//          ftempreco_ncfit_post->Scale(nc_shift);
//          ftemptrue->Reset();
//          //Oscillate true cc's selected as nc's after efficiency and turn them to reco
//          //oscweight(fnooscar_ccselnc_post,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
//          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnc_post,ftemptrue,14,NCType::kCC);
//          truetoreco(ftemptrue,*ftemprecocc_ncfit_post,matarrayb_post[ii]);
//          ftemptrue->Reset();

//          // TAUS  HERE
//          if(fDoMC==0){
//            //Oscillate true cc's selected as cc's after efficiency and turn them to reco
//            // LAST 11 IS FOR CC TAU
//            //oscweight(fnooscar_ccselcct_post,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
//            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcct_post,ftemptrue,16,NCType::kCC);
//            truetoreco(ftemptrue,*ftemprecotau_post,matarrayt_post[ii]);
//            ftemptrue->Reset();

//            //Oscillate true cc's selected as nc's after efficiency and turn them to reco
//            //oscweight(fnooscar_ccselnct_post,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
//            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnct_post,ftemptrue,16,NCType::kCC);
//            truetoreco(ftemptrue,*ftemprecotau_ncfit_post,matarraybt_post[ii]);
//            ftemptrue->Reset();
//          }

            chisqmin_norm = 1.e12;
            // NORMALIZATION // kk
            for(int kk = nr0; kk < nr1; kk++){
              norm = lnr + kk*size_nr;
              // CC
              ftempreconc->Scale(norm);
              ftempreco->Scale(norm);
//            ftempreconc_post->Scale(norm);
//            ftempreco_post->Scale(norm);
              if(fDoMC==0){
                ftemprecotau->Scale(norm);
//              ftemprecotau_post->Scale(norm);
              }
              // NC
              ftemprecocc_ncfit->Scale(norm);
              ftempreco_ncfit->Scale(norm);
//            ftemprecocc_ncfit_post->Scale(norm);
//            ftempreco_ncfit_post->Scale(norm);
              if(fDoMC==0){
                ftemprecotau_ncfit->Scale(norm);
//              ftemprecotau_ncfit_post->Scale(norm);
              }
              //ADD THEM ALL CC
              ftempreco->Add(ftempreconc);
//            ftempreco_post->Add(ftempreconcp);
              //ADD THEM ALL NC
              ftempreco_ncfit->Add(ftemprecocc_ncfit);
//            ftempreco_ncfit_post->Add(ftemprecocc_ncfitp);
              if(fDoMC == 0){
                ftempreco->Add(ftemprecotau);
//              ftempreco_post->Add(ftemprecotaup);
                ftempreco_ncfit->Add(ftemprecotau_ncfit);
//              ftempreco_ncfit_post->Add(ftemprecotau_ncfitp);
              }

              //REBIN
              Rebin(ftempreco, *ftempreco_rebinned);
              Rebin(ftempreconc, *ftempreconc_rebinned);
//            Rebin(ftempreco_post, *ftempreco_rebinned_post);
//            Rebin(ftempreconc_post, *ftempreconc_rebinned_post);
              Rebin(ftempreco_ncfit, *ftempreco_ncfit_rebinned);
              Rebin(ftemprecocc_ncfit, *ftemprecocc_ncfit_rebinned);
//            Rebin(ftempreco_ncfit_post, *ftempreco_ncfit_rebinned_post);
//            Rebin(ftemprecocc_ncfit_post, *ftemprecocc_ncfit_rebinned_post);

              // CHISQUARE
              chisqtemp1 = chisq(foscreco_rebinned, ftempreco_rebinned, ftempreconc_rebinned);
              chisqtemp1 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
                              +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
                              +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
//            chisqtemp2 = chisq(foscreco_rebinned_post, ftempreco_rebinned_post, ftempreconc_rebinned_post);
//            chisqtemp2 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
//                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
//                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
              chisqtemp3 = chisq(foscreco_ncfit_rebinned, ftempreco_ncfit_rebinned, ftemprecocc_ncfit_rebinned);
//            chisqtemp3 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
//                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
//                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
//            chisqtemp4 = chisq(foscreco_ncfit_rebinned_post, ftempreco_ncfit_rebinned_post, ftemprecocc_ncfit_rebinned_post);
//            chisqtemp4 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
//                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
//                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );

              chisqtemp1 += chisqtemp3;//chisqtemp2 + chisqtemp3 + chisqtemp4;

              if(chisqtemp1 < chisqmin_norm) chisqmin_norm = chisqtemp1;

              if(chisqtemp1 < chisqmin){
                chisqmin = chisqtemp1;
                dmsqmin  = dmsqfit;
                s2tmin   = s2tfit;
                s2smin   = s2sfit;
                norm_min = norm;
                nc_min   = nc_shift;
                sh_min   = enu_shift;
                best_enu = ii;
                best_nc  = jj;
              }
              // Un-Normalize
              if(fDoMC == 0){
                ftempreco->Add(ftemprecotau,-1);
                ftemprecotau->Scale(1./norm);
//              ftempreco_post->Add(ftemprecotau_post,-1);
//              ftemprecotau_post->Scale(1./norm);
                ftempreco_ncfit->Add(ftemprecotau_ncfit,-1);
                ftemprecotau_ncfit->Scale(1./norm);
//              ftempreco_ncfit_post->Add(ftemprecotau_ncfit_post,-1);
//              ftemprecotau_ncfit_post->Scale(1./norm);
              }
              ftempreco->Add(ftempreconc,-1);
              ftempreco->Scale(1./norm);
              ftempreconc->Scale(1./norm);
//            ftempreco_post->Add(ftempreconc_post,-1);
//            ftempreco_post->Scale(1./norm);
//            ftempreconc_post->Scale(1./norm);
              ftempreco_ncfit->Add(ftemprecocc_ncfit,-1);
              ftempreco_ncfit->Scale(1./norm);
              ftemprecocc_ncfit->Scale(1./norm);
//            ftempreco_ncfit_post->Add(ftemprecocc_ncfit_post,-1);
//            ftempreco_ncfit_post->Scale(1./norm);
//            ftemprecocc_ncfit_post->Scale(1./norm);

            }// NORMALIZATION

            chisqarraymins[i][k][j] = chisqmin_norm;

            MSG("NCExtrapolationNS", Msg::kInfo) << " (dm^2, sin^2 2theta, fsterile) = ("
                                                 << dmsqfit << "," << s2tfit << ","
                                                 << s2sfit << ") chi^2 = "
                                                 << chisqmin_norm << endl;
          }// fs
        }// sin2theta
      }// dm2

      //Unnormalized NC
      ftempreconc->Scale(1./nc_shift);
      ftempreco_ncfit->Scale(1./nc_shift);
//       ftempreconc_post->Scale(1./nc_shift);
//       ftempreco_ncfit_post->Scale(1./nc_shift);

    }// NC LOOP
  }// SHOWER LOOP


  //set the best fit values
  BestUMu3Sqr() = s2tmin;
  BestUS3Sqr() = s2smin;
  BestDeltaMSqr() = dmsqmin;
  fMinChiSqr = chisqmin;

  cout <<"Bestfit values stored as Oscillation Parameters!: "
       <<"\nUMu3SqrTheta = " << BestUMu3Sqr()
       <<"\nUS3Sqr = " << BestUS3Sqr()
       <<"\nDeltaMSqr = " << BestDeltaMSqr()
       <<"\nMinChiSqr = " << fMinChiSqr
       << endl;

//   fBeams[farBeamI].SetOscillationParameters(fBestUMu3Sqr, fBestDeltaMSqr, fBestUS3Sqr);
//   fBeams[farBeamII].SetOscillationParameters(fBestUMu3Sqr, fBestDeltaMSqr, fBestUS3Sqr);

  NCExtrapolation::FillDeltaChiSqrMaps(chisqarraymins);

  // Now fill chisq
  for(int i = 0; i < maxdm; i++){// dm2
    dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;

    for(int j = 0; j < maxsin; j++){ //sin2theta
      s2tfit = (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;

      double chismin =1e9;
      int    kmin    =-1;
      for(int k = 0; k < maxfs; k++){ //fsterile
        s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;

        if(chisqarraymins[i][k][j]<chismin){
          chismin = chisqarraymins[i][k][j];
          kmin = k;
        }
      } // fs
      fchisqall->Fill(s2tfit, dmsqfit, chismin);
    } // sin2theta
  } // dm2

  double chismin = 1e9;
  int    jmin    = -1;

  // Now fill chisq other projection
  for(int i = 0; i < maxdm; i++){// dm2
    dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;

    for(int k = 0; k < maxfs; k++){ //fsterile
      s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;

      chismin = 1e9;
      jmin    = -1;

      for(int j = 0; j < maxsin; j++){ //sin2theta
        s2tfit = (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;

        if(chisqarraymins[i][k][j] < chismin){
          chismin = chisqarraymins[i][k][j];
          jmin = j;
        }
      } // sin2theta
      fchisqallnc->Fill(s2sfit, dmsqfit, chismin);
    } // fs
  } // dm2

  TString name = "";
  TString junk = "";
  TString domch = "data";
  TString dodecayh = "osc";

  if(fDoMC==1) domch = "mc";
  if(fDecay==1) dodecayh = "dec";

  int sele = n1d/2;
  int selnc = ncd/2;
  enu_shiftmin = 1.;
  nc_shiftmin  = 1.;
  if(fSystematicsToUse[NCType::kShowerEnergy]){
    sele = best_enu;
    enu_shiftmin = les+best_enu*size_e;
  }
  if(fSystematicsToUse[NCType::kNCBackground]){
    selnc = best_nc;
    nc_shiftmin = lnc+best_nc *size_nc;
  }

  if(fDoMC==1){
    name = "res_fit_";
    name += "_" + junk + "_" + domch + "_" + dodecayh;
    name += "new.txt";
  }
  if(fDoMC==0){
    name = "res_fit_" + domch + "_"  + dodecayh;
    name += "_test_post_part_new.txt";
  }

  //Temp! - removal of systematic loops, setting normalizations to 1
  norm_min = 1.0;
  enu_shiftmin = 1.0;
  nc_shiftmin = 1.0;

  ofstream text_file(name);
  text_file << " Best fit at dmsq= " << dmsqmin
            << " sin2theta= "       << s2tmin
            << " fs " << s2smin
            << "  chisq=" << chisqmin/ndf
            << " norm_min " << norm_min
            << " enu_shiftmin " << enu_shiftmin
            << " nc_shiftmin " << nc_shiftmin << endl;
  text_file.close();

  MSG("NCExtrapolationNS", Msg::kInfo) << " Best fit at dmsq = " << dmsqmin
                                       << " sin2theta = " << s2tmin
                                       << " fs = " << s2smin << "  chisq = "
                                       << chisqmin << " ndf = "  << ndf << " norm_min "
                                       << norm_min << " enu_shiftmin "
                                       << enu_shiftmin << "nc_shiftmin " << nc_shiftmin
                                       <<  endl;
  MSG("NCExtrapolationNS", Msg::kInfo) << " Best enu integer= " << best_enu
                                       << " best nc integer "  << best_nc << endl;

  // Make best fit histo//
  // For that we need best shower energy, best nc value,
  //best normalization value and the best fit points//
  fnooscar->Reset();
//   fnooscar_post->Reset();
  fnooscar_ccselcc->Reset();
//   fnooscar_ccselcc_post->Reset();
  fnooscar_ccselnc->Reset();
//   fnooscar_ccselnc_post->Reset();
  fnooscar_ncselnc->Reset();
//   fnooscar_ncselnc_post->Reset();
  fnooscar_ncselcc->Reset();
//   fnooscar_ncselcc_post->Reset();
  // taus
  fnooscar_ccselcct->Reset();
//   fnooscar_ccselcct_post->Reset();
  fnooscar_ccselnct->Reset();
//   fnooscar_ccselnct_post->Reset();

  //Copying noosctruearray for best showerenergy and nc values
//Temp disabling to check fitting, no looping over showerenergy and nc values at the moment   noosctruearray[sele][selnc]->Copy(*fnooscar);
  fallccmcf->Copy(*fnooscar); //temp using true FD spectrum (before pid efficiency corr) instead of prediction!
//   noosctruearraypost[sele][selnc]->Copy(*fnooscar_post);

  // PRE
  fnooscar_ccselcc->Multiply(fnooscar,feffic);
  fnooscar_ccselnc->Multiply(fnooscar,fefficbg);
  fnooscar_ncselnc->Multiply(fnooscar,fratccnc);
  fnooscar_ncselcc->Multiply(fnooscar,fratccnc);
  fnooscar_ncselnc->Multiply(fefficnc);
  fnooscar_ncselcc->Multiply(fefficncbg);

  // POST
//   fnooscar_ccselcc_post->Multiply(fnooscar_post,fefficp);
//   fnooscar_ccselnc_post->Multiply(fnooscar_post,fefficbgp);
//   fnooscar_ncselnc_post->Multiply(fnooscar_post,fratccncp);
//   fnooscar_ncselcc_post->Multiply(fnooscar_post,fratccncp);
//   fnooscar_ncselnc_post->Multiply(fefficncp);
//   fnooscar_ncselcc_post->Multiply(fefficncbgp);

  // taus
  if(fDoMC==0){
    // PRE
    fnooscar_ccselcct->Multiply(fnooscar,frattau);
    fnooscar_ccselcct->Multiply(feffict);
    fnooscar_ccselnct->Multiply(fnooscar,frattau);
    fnooscar_ccselnct ->Multiply(fefficbgt);

    // POST
//     fnooscar_ccselcct_post->Multiply(fnooscar_post,frattaup);
//     fnooscar_ccselcct_post->Multiply(feffictp);
//     fnooscar_ccselnct_post->Multiply(fnooscar_post,frattaup);
//     fnooscar_ccselnct_post->Multiply(fefficbgtp);

  }

  //PRE SHUTDOWN

  for(int dmb = 0; dmb < fnbins; ++dmb){
    sinSqrDeltaMSqr.push_back(TMath::Power(TMath::Sin(NCType::k127*dmsqfit
                                                      *NCType::kBaseLineFar/fnooscreco->GetBinCenter(dmb+1)),2.));
  }//end loop over energy bins

  //Oscillate true cc's selected as cc's after efficiency and turn them to reco
  //oscweight(fnooscar_ccselcc,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
  oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcc,ftemptrue,14,NCType::kCC);
  truetoreco(ftemptrue,*fbestfitreco,matarray[sele]);
  ftemptrue->Reset();
  //Oscillate true nc's selected as cc's after efficiency and turn them to reco
  //oscweight(fnooscar_ncselcc,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
  oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselcc,ftemptrue,14,NCType::kNC);
  truetoreco(ftemptrue,*fbestfitreconc,matarrayncb[sele]);
  fbestfitreconc->Scale(nc_shiftmin);
  ftemptrue->Reset();
  //Oscillate true nc's selected as nc's after efficiency and turn them to reco
  //oscweight(fnooscar_ncselnc,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
  oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselnc,ftemptrue,14,NCType::kNC);
  truetoreco(ftemptrue,*fbestfitreco_ncfit,matarraync[sele]);
  fbestfitreco_ncfit->Scale(nc_shiftmin);
  ftemptrue->Reset();
  //Oscillate true cc's selected as nc's after efficiency and turn them to reco
  //oscweight(fnooscar_ccselnc,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
  oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnc,ftemptrue,14,NCType::kCC);
  truetoreco(ftemptrue,*fbestfitrecocc_ncfit,matarrayb[sele]);
  ftemptrue->Reset();

  // TAUS  HERE
  if(fDoMC==0){
    //Oscillate true cc's selected as cc's after efficiency and turn them to reco
    // LAST 11 IS FOR CC TAU
    //oscweight(fnooscar_ccselcct,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
    oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcct,ftemptrue,16,NCType::kCC);
    truetoreco(ftemptrue,*fbestfitrecotau,matarrayt[sele]);
    ftemptrue->Reset();

    //Oscillate true cc's selected as nc's after efficiency and turn them to reco
    //oscweight(fnooscar_ccselnct,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
    oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnct,ftemptrue,16,NCType::kCC);
    truetoreco(ftemptrue,*fbestfitrecotau_ncfit,matarraybt[sele]);
    ftemptrue->Reset();
  }

//   //POST SHUTDOWN
//   //Oscillate true cc's selected as cc's after efficiency and turn them to reco
//   //oscweight(fnooscar_ccselcc_post,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
//   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcc_post,ftemptrue,14,NCType::kCC);
//   truetoreco(ftemptrue,*fbestfitreco_post,matarray_post[sele]);
//   ftemptrue->Reset();
//   //Oscillate true nc's selected as cc's after efficiency and turn them to reco
//   //oscweight(fnooscar_ncselcc_post,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
//   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselcc_post,ftemptrue,14,NCType::kNC);
//   truetoreco(ftemptrue,*fbestfitreconc_post,matarrayncb_post[sele]);
//   fbestfitreconc_post->Scale(nc_shiftmin);
//   ftemptrue->Reset();
//   //Oscillate true nc's selected as nc's after efficiency and turn them to reco
//   //oscweight(fnooscar_ncselnc_post,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
//   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselnc_post,ftemptrue,14,NCType::kNC);
//   truetoreco(ftemptrue,*fbestfitreco_ncfit_post,matarraync_post[sele]);
//   fbestfitreco_ncfit_post->Scale(nc_shiftmin);
//   ftemptrue->Reset();
//   //Oscillate true cc's selected as nc's after efficiency and turn them to reco
//   //oscweight(fnooscar_ccselnc_post,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
//   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnc_post,ftemptrue,14,NCType::kCC);
//   truetoreco(ftemptrue,*fbestfitrecocc_ncfit_post,matarrayb_post[sele]);
//   ftemptrue->Reset();

//   // TAUS  HERE
//   if(fDoMC==0){
//     //Oscillate true cc's selected as cc's after efficiency and turn them to reco
//     //oscweight(fnooscar_ccselcct_post,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
//     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcct_post,ftemptrue,16,NCType::kCC);
//     truetoreco(ftemptrue,*fbestfitrecotau_post,matarrayt_post[sele]);
//     ftemptrue->Reset();

//     //Oscillate true cc's selected as nc's after efficiency and turn them to reco
//     //oscweight(fnooscar_ccselnct_post,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
//     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnct_post,ftemptrue,16,NCType::kCC);
//     truetoreco(ftemptrue,*fbestfitrecotau_ncfit_post,matarraybt_post[sele]);
//     ftemptrue->Reset();
//   }

  // CC
  fbestfitreconc->Scale(norm_min);
  fbestfitreco->Scale(norm_min);
//   fbestfitreconc_post->Scale(norm_min);
//   fbestfitreco_post->Scale(norm_min);
  fbestfitrecocc_ncfit->Scale(norm_min);
  fbestfitreco_ncfit->Scale(norm_min);
  fbestfitrecocc_ncfit->Scale(norm_min);
//   fbestfitreco_ncfit_post->Scale(norm_min);
  fbestfitreco->Add(fbestfitreconc);
//   fbestfitreco_post->Add(fbestfitreconcp);
  fbestfitreco_ncfit->Add(fbestfitrecocc_ncfit);
//   fbestfitreco_ncfit_post->Add(fbestfitrecocc_ncfitp);
  if(fDoMC==0){
    fbestfitrecotau->Scale(norm_min);
//     fbestfitrecotau_post->Scale(norm_min);
    fbestfitrecotau_ncfit->Scale(norm_min);
//     fbestfitrecotau_ncfit_post->Scale(norm_min);
    fbestfitreco->Add(fbestfitrecotau);
//     fbestfitreco_post->Add(fbestfitrecotaup);
    fbestfitreco_ncfit->Add(fbestfitrecotau_ncfit);
//     fbestfitreco_ncfit_post->Add(fbestfitrecotau_ncfitp);
  }

  //REBIN
  Rebin(fbestfitreco, *fbestfitreco_rebinned);
  Rebin(fbestfitreconc, *fbestfitreconc_rebinned);
//   Rebin(fbestfitreco_post, *fbestfitreco_rebinned_post);
//   Rebin(fbestfitreconc_post, *fbestfitreconc_rebinned_post);
  Rebin(fbestfitreco_ncfit, *fbestfitreco_ncfit_rebinned);
  Rebin(fbestfitrecocc_ncfit, *fbestfitrecocc_ncfit_rebinned);
//   Rebin(fbestfitreco_ncfit_post, *fbestfitreco_ncfit_rebinned_post);
//   Rebin(fbestfitrecocc_ncfit_post, *fbestfitrecocc_ncfit_rebinned_post);
  if(fDoMC==0){
    Rebin(fbestfitrecotau, *fbestfitrecotau_rebinned);
    Rebin(fbestfitrecotau_ncfit, *fbestfitrecotau_ncfit_rebinned);
//     Rebin(fbestfitrecotau_post, *fbestfitrecotau_rebinned_post);
//     Rebin(fbestfitrecotau_ncfit_post, *fbestfitrecotau_ncfit_rebinned_post);
  }

  // PREDICTION UNOSCILLATED //
  fnooscar->Copy(*fnoosctrue);
//   fnooscar_post->Copy(*fnoosctruep);
  fnoosctrue_ncfit->Multiply(fnooscar,fratccnc);
//   fnoosctrue_ncfit_post->Multiply(fnooscar_post,fratccncp);

  //PRE SHUTDOWN
  truetoreco(fnooscar_ccselcc,*fnooscreco,  matarray[sele]);
  truetoreco(fnooscar_ncselcc,*fnooscreconc,matarrayncb[sele]);
  fnooscreconc->Scale(nc_shiftmin);
  truetoreco(fnooscar_ncselnc,*fnooscreco_ncfit,matarraync[sele]);
  fnooscreconc->Scale(nc_shiftmin);
  truetoreco(fnooscar_ccselnc,*fnooscrecocc_ncfit, matarrayb[sele]);

  //POST SHUTDOWN
  truetoreco(fnooscar_ccselcc_post,*fnooscreco_post,matarray_post[sele]);
  truetoreco(fnooscar_ncselcc_post,*fnooscreconc_post,matarrayncb_post[sele]);
//   fnooscreconc_post->Scale(nc_shiftmin);
  truetoreco(fnooscar_ncselnc_post,*fnooscreco_ncfit_post,matarraync_post[sele]);
//   fnooscreconc_post->Scale(nc_shiftmin);
  truetoreco(fnooscar_ccselnc_post,*fnooscrecocc_ncfit_post,matarrayb_post[sele]);

  // CC
  fnooscreconc->Scale(norm_min);
  fnooscreco->Scale(norm_min);
//   fnooscreconc_post->Scale(norm_min);
  fnooscreco_post->Scale(norm_min);
  fnooscrecocc_ncfit->Scale(norm_min);
  fnooscreco_ncfit->Scale(norm_min);
//   fnooscrecocc_ncfit_post->Scale(norm_min);
//   fnooscreco_ncfit_post->Scale(norm_min);
  fnooscreco->Add(fnooscreconc);
//   fnooscreco_post->Add(fnooscreconcp);
  fnooscreco_ncfit->Add(fnooscrecocc_ncfit);
//   fnooscreconc_post->Add(fnooscrecocc_ncfitp);
  if(fDoMC==0){
    fnooscrecotau->Scale(norm_min);
//     fnooscrecotau_post->Scale(norm_min);
    fnooscrecotau_ncfit->Scale(norm_min);
//     fnooscrecotau_ncfit_post->Scale(norm_min);
    fnooscreco->Add(fnooscrecotau);
//     fnooscreco_post->Add(fnooscrecotaup);
    fnooscreco_ncfit->Add(fnooscrecotau_ncfit);
//     fnooscreconc_post->Add(fnooscrecotau_ncfitp);
  }

  //REBIN
  Rebin(fnooscreco,          *fnooscreco_rebinned);
  Rebin(fnooscreconc,        *fnooscreconc_rebinned);
//   Rebin(fnooscreco_post,         *fnooscreco_rebinned_post);
//   Rebin(fnooscreconc_post,       *fnooscreconc_rebinned_post);
  Rebin(fnooscreco_ncfit,       *fnooscreco_ncfit_rebinned);
  Rebin(fnooscrecocc_ncfit,        *fnooscrecocc_ncfit_rebinned);
//   Rebin(fnooscreco_ncfit_post,      *fnooscreco_ncfit_rebinned_post);
//   Rebin(fnooscrecocc_ncfit_post,       *fnooscrecocc_ncfit_rebinned_post);
  if(fDoMC==0){
    Rebin(fnooscrecotau, *fnooscrecotau_rebinned);
    Rebin(fnooscrecotau_ncfit, *fnooscrecotau_ncfit_rebinned);
//     Rebin(fnooscrecotau_post, *fnooscrecotau_rebinned_post);
//     Rebin(fnooscrecotau_ncfit_post, *fnooscrecotau_ncfit_rebinned_post);
  }

//   //clear all the far beam MC histograms and then add the
//   //results from the extrapolations to the appropriate ones
// LLH uncommented
  int farBeams[] = {farBeamI, farBeamII};
  for(int b = 0; b < 2; ++b){
     for(int i = 0; i < 2; ++i){
        fBeams[farBeams[b]].GetMCFitHistogram(i)->Reset();
        fBeams[farBeams[b]].GetMCBackgroundHistogram(i)->Reset();
        fBeams[farBeams[b]].GetMCFitNuMuToNuTauHistogram(i)->Reset();
        fBeams[farBeams[b]].GetMCHistogram(i)->Reset();
        fBeams[farBeams[b]].GetMCBackgroundHistogram(i)->Reset();
        fBeams[farBeams[b]].GetMCNoFitNuMuToNuTauHistogram(i)->Reset();
     }
  }

  float scaleFactor = 1.; //what is this?  LLH

//LLH uncommented
   //put the fit spectra into the beam objects
   fBeams[farBeamI].GetMCFitHistogram(NCType::kCC)->Add(fbestfitreco_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCFitHistogram(NCType::kNC)->Add(fbestfitreco_ncfit_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kCC)->Add(fbestfitreconc_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kNC)->Add(fbestfitrecocc_ncfit_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCFitNuMuToNuTauHistogram(NCType::kCC)->Add(fbestfitrecotau_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCFitNuMuToNuTauHistogram(NCType::kNC)->Add(fbestfitrecotau_ncfit_rebinned, 1./scaleFactor);

   //and the nooscillation histograms
   fBeams[farBeamI].GetMCHistogram(NCType::kCC)->Add(fnooscreco_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCHistogram(NCType::kNC)->Add(fnooscreco_ncfit_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kCC)->Add(fnooscreconc_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kNC)->Add(fnooscrecocc_ncfit_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCNoFitNuMuToNuTauHistogram(NCType::kCC)->Add(fnooscrecotau_rebinned, 1./scaleFactor);
   fBeams[farBeamI].GetMCNoFitNuMuToNuTauHistogram(NCType::kNC)->Add(fnooscrecotau_ncfit_rebinned, 1./scaleFactor);

   //put the fit spectra into the beam objects
   fBeams[farBeamII].GetMCFitHistogram(NCType::kCC)->Add(fbestfitreco_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCFitHistogram(NCType::kNC)->Add(fbestfitreco_ncfit_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kCC)->Add(fbestfitreconc_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kNC)->Add(fbestfitrecocc_ncfit_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCFitNuMuToNuTauHistogram(NCType::kCC)->Add(fbestfitrecotau_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCFitNuMuToNuTauHistogram(NCType::kNC)->Add(fbestfitrecotau_ncfit_rebinned_post, 1./scaleFactor);

   //and the nooscillation histograms
   fBeams[farBeamII].GetMCHistogram(NCType::kCC)->Add(fnooscreco_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCHistogram(NCType::kNC)->Add(fnooscreco_ncfit_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kCC)->Add(fnooscreconc_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kNC)->Add(fnooscrecocc_ncfit_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCNoFitNuMuToNuTauHistogram(NCType::kCC)->Add(fnooscrecotau_rebinned_post, 1./scaleFactor);
   fBeams[farBeamII].GetMCNoFitNuMuToNuTauHistogram(NCType::kNC)->Add(fnooscrecotau_ncfit_rebinned_post, 1./scaleFactor);

  return;
}

//----------------------------------------------------------------------
//write out histograms etc used in the fit
void NCExtrapolationNS::WriteResources()
{

  //fill the histograms in the NCBeam objects in fBeams vector
  //before calling WriteResources in the base class

   cout <<"Writing!!!!!!" << endl;

  NCExtrapolation::WriteResources();

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* fileDir = gDirectory;
  TDirectory* saveDir = gDirectory;
  fileDir->cd("plots"+NCType::kExtrapolationNames[fExtrapolationType]);

  fTrueND->Write();
  fTrueNDAfterRT->Write();
  fTrueNDAfterRT_noeffcorr->Write(); //LLH adding
  fRecoND->Write();
  fTrueFD->Write();
  fTrueFD_nc->Write();
  fRecoFD->Write();
  fRecoFD_nc->Write();
  fTrueFDAfterBM->Write();
  fTrueFDAfterBM_noacccorr->Write(); //LLH addng
  fTrueFDAfterBM_posteffcorr->Write(); //LLH addng
  fTrueFDAfterOsc->Write(); //LLH addng
  fRecoFDAfterBM->Write();

  //LLH adding
  frecotruemat->Write();
  fefficiencyn->Write();
  fdatselND->Write();
  fdatselFD->Write();
  fdatselFD_nc->Write();

  //true to reco matricies
  for(uint i = 0; i < matarray.size(); ++i){
    matarray[i]->Write();
    matarraync[i]->Write();
    matarrayb[i]->Write();
    matarrayncb[i]->Write();
    //TAUS
    matarrayt[i]->Write();
    matarraybt[i]->Write();

    matarray_post[i]->Write();
    matarraync_post[i]->Write();
    matarrayb_post[i]->Write();
    matarrayncb_post[i]->Write();
    //TAUS
    matarrayt_post[i]->Write();
    matarraybt_post[i]->Write();

    for(uint j = 0; j < noosctruearray[i].size(); ++j){
      noosctruearray[i][j]->Write();
      noosctruearraypost[i][j]->Write();
    }
  }//end writing true to reco matricies


  fnoosctrue->Write();
  fnooscreco->Write();
  fbestfitreco->Write();
  fbestfitreconc->Write();
  fbestfitrecotau->Write();

  //LLH adding
  fnooscar->Write();
  fnooscar_ccselcc->Write();
  fnooscar_ccselnc->Write();
  fnooscar_ncselnc->Write();
  fnooscar_ncselcc->Write();

  fnooscar_ccselcct->Write();
  fnooscar_ccselnct->Write();

  foscreco_rebinned->Write();
  fnooscreco_rebinned->Write();
  fbestfitreco_rebinned->Write();
  fbestfitreconc_rebinned->Write();
  fbestfitrecotau_rebinned->Write();
  //
  fnoosctrue_ncfit->Write();
  fnooscreco_ncfit->Write();
  fnooscrecocc_ncfit->Write();
  fnooscrecotau_ncfit->Write();
  fbestfitreco_ncfit->Write();
  fbestfitrecocc_ncfit->Write();
  fbestfitrecotau_ncfit->Write();

  foscreco_ncfit_rebinned->Write();

  fnooscreco_ncfit_rebinned->Write();
  fbestfitreco_ncfit_rebinned->Write();
  fbestfitrecocc_ncfit_rebinned->Write();
  fbestfitrecotau_ncfit_rebinned->Write();

  // POst

  fnoosctrue_post->Write();
  fnooscreco_post->Write();
  fbestfitreco_post->Write();
  fbestfitreconc_post->Write();
  fbestfitrecotau_post->Write();

  foscreco_rebinned_post->Write();

  fnooscreco_rebinned_post->Write();
  fbestfitreco_rebinned_post->Write();
  fbestfitreconc_rebinned_post->Write();
  fbestfitrecotau_rebinned_post->Write();

  fnoosctrue_ncfit_post->Write();
  fnooscreconc_post->Write();
  fbestfitreco_ncfit_post->Write();
  fbestfitrecocc_ncfit_post->Write();
  fbestfitrecotau_ncfit_post->Write();

  foscreco_ncfit_rebinned_post->Write();

  fnooscreco_ncfit_rebinned_post->Write();
  fbestfitreco_ncfit_rebinned_post->Write();
  fbestfitrecocc_ncfit_rebinned_post->Write();
  fbestfitrecotau_ncfit_rebinned_post->Write();

  //
  fchisqall->Write();
  fchisqallnc->Write();
  fratccnc->Write();

  fcorf->Write();
  fcor->Write();
  fcorfall->Write();
  fcorall->Write();
  fSingleMCNearNuCC->Write();
  fSingleMCNearAllCC->Write();
  fSingleMCFarNuCC->Write();
  fSingleMCFarAllCC->Write();
  fMCNearNuCC->Write();
  fMCNearAllCC->Write();
  fMCFarNuCC->Write();
  fMCFarAllCC->Write();
  fallccmc->Write();
  fallccmcf->Write();
  fallncmcf->Write();

  fpurityn->Write();
  feffic->Write();
  fefficnc->Write();
  fefficbg->Write();
  fefficncbg->Write();
  //fratccnc->Write(); //LLH oops in here twice!
  frattau->Write();
  feffic_numcheck->Write(); //LLH adding

  saveDir->cd();

  return;
}

//***************************************************************************************
void NCExtrapolationNS::MakePrediction(TH1F *datsel,
                                       TH1F* purityn,
                                       TH1F* efficiencyn,
                                       TH2F *matn_norm,
                                       TH1F* prediction,
                                       int runType){

  MSG("NCExtrapolationNS", Msg::kInfo) << " in MakePrediction" << endl;

  prediction->Reset();

  TH2D *bmat_norm_r2 = fbmat_norm_r2;
  if(runType == 1) bmat_norm_r2 = fbmat_norm_r2_post; //LLH should this be runType == NCType::kRunII?? - fixme

  double ntrue[fnbins] = {0.};
  double ftrue[fnbins] = {0.};
  double eff[fnbins] = {0.};

  // FIDUCIAL MASS DEFINITIONS
  double scalef = kscalef;
  double mass_f = kmass_f;
  double mass_n = kmass_n;

//this was old fv for PRL (doMC = 0 is data, doMC=1 is mc)
//   if(fDoMC==-1){ // !!!!!!!!! here change
//     mass_f    = 5.23*1e3+53*484*0.01*1.032; // tons
//     mass_n    = 3.141592654*(1.)*(1.)*((84-17)*0.0254*7.874+(84-17)*0.01*1.032);
//     scalef    = 1.00;
//   }
//commented out, this is the old cc fv volume - LLH

  double  scalemass = mass_f/mass_n;
  cout <<"scalemass = " << scalemass
       <<"\nmass_f = " << mass_f
       <<"\nmass_n = " << mass_n
       << endl;

  TH1F *datsel_cor = new TH1F("datsel_cor","",fnbins,fbinvect);
  // NEAR TRUE SPECTRA
  TH1F *allccmc_obt_2d = new TH1F("allccmc_obt_2d","",fnbins,fbinvect);

  // STEP 1 : CORRECT FOR PURITY THE "DATA SELECTED SPECTGRUM"
  datsel_cor->Multiply(datsel,purityn);

  // STEP 2 : TRANSFORM RECO = > TRUE & CORRECT FOR EFFICIENCY
  float k = 0.;
  for(int i = 0; i < fnbins; i++){
    ntrue[i] = 0.;
    eff[i] = efficiencyn->GetBinContent(i+1);
    k = datsel_cor->GetBinCenter(i+1);

    // TRANSFORMATION USING 2D MATRIX
    for(int j = 0; j < fnbins; j++){
      ntrue[i] += matn_norm->GetBinContent(j+1,i+1)*datsel_cor->GetBinContent(j+1);
    }
    if(eff[i] > 0.){
       fTrueNDAfterRT_noeffcorr->Fill(k, ntrue[i]); //LLH adding
       ntrue[i] /= eff[i];   //efficiency correction!
       allccmc_obt_2d->Fill(k, ntrue[i]);
       fTrueNDAfterRT->Fill(k, ntrue[i]);
    }

//LLH commenting out
//     MAXMSG("NCExtrapolationNS", Msg::kInfo, fnbins)
//       << "near reco to true: bin center = " << k
//       << ", datsel = " << datsel->GetBinContent(i+1)
//       << ", datsel_cor (purity corr) = " << datsel_cor->GetBinContent(i+1)
//       << ", ntrue (eff and reco->true corr) = " << ntrue[i] << endl;

  }//end loop over bins

  // STEP 3 EXTRAPOLATE & CORRECT FOR ACCEPTANCE
  double cor_factf = 0.;

  for(int i = 0; i < fnbins; i++){
    // Correction factor for acceptance
    cor_factf = fcorf->GetBinContent(i+1);
    ftrue[i] = 0.;

    // TRANSFORMATION USING 2D GNUMI MATRIX
    for(int j = 0; j < fnbins; j++) {
      ftrue[i] += bmat_norm_r2->GetBinContent(j+1,i+1)*allccmc_obt_2d->GetBinContent(j+1);
    } //LLH added braces for readability

    //LLH adding
    k = fcorf->GetBinCenter(i+1);
    fTrueFDAfterBM_noacccorr->Fill(k, ftrue[i]);
    //end LLH added

//LLH    cout <<"cor_factf = " << cor_factf << endl;

    prediction->SetBinContent(i+1,ftrue[i]*cor_factf);  //applying correction for FD acceptance and event recon.
  }

  prediction->Scale(scalemass*scalef);

  datsel_cor->Delete();
  allccmc_obt_2d->Delete();
}

//----------------------------------------------------------------------
void NCExtrapolationNS::RecoTruePur(TH2F*recotruemat1,
                                    TH1F *purityn1,
                                    TH1F *efficiencyn1,
                                    double enu_shift,
                                    int beam)
{
   cout <<"In RecoTruePur! " << endl;

  recotruemat1->Reset();
  purityn1->Reset();
  efficiencyn1->Reset();

  TH1F selcc_true("selcc_true","",fnbins,fbinvect);
  TH1F selcc_reco("selcc_reco","",fnbins,fbinvect);
  TH1F sel_reco("sel_reco","",fnbins,fbinvect);

  double mc_weif=1.;
  double shower = 0.;
  double mumom = 0.;
  double reco_ene = 0.;
  int    flavour = 0;
  double true_enu = 0.;
  //double annpid = 0.;
  //int cc_nc = 0;
  //int initial_state = 0;

  //size of the MC arrays
  int sizeSignal =  0;
  int sizeBG = 0;

  //loop over bins in the beam
  //only use CC events in this stage
  NCEnergyBin *bin = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
    sizeSignal = bin->GetMCSignalVectorSize();
    sizeBG = bin->GetMCBackgroundVectorSize();
    for(int j = 0; j < sizeSignal; ++j){
      bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      //only use nu_mu, no nu_mubar
      if(flavour == 14){
        selcc_true.Fill(true_enu,mc_weif);
        recotruemat1->Fill(reco_ene,true_enu,mc_weif);
        selcc_reco.Fill(reco_ene,mc_weif);
        fTrueND->Fill(true_enu, mc_weif); } //LLH temp to check datsel against fRecoND
        fRecoND->Fill(reco_ene, mc_weif);
//LLH temp default brace location     }
        sel_reco.Fill(reco_ene,mc_weif); //LLH sel_reco includes numu and numubar
    }//end loop over signal events

    for(int j = 0; j < sizeBG; ++j){
      bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      sel_reco.Fill(reco_ene,mc_weif);
      fRecoND->Fill(reco_ene, mc_weif); //LLH temp to check datsel against fRecoND
    }//end loop over background events

    //LLH adding
    int size_elec = bin->GetMCNuEToNuEVectorSize();
    int size_tau = bin->GetMCNuMuToNuTauVectorSize(); //there should be no taus in ND!
    for(int j = 0; j < size_tau; ++j){
       bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       shower *= enu_shift;
       reco_ene = shower + mumom;
       sel_reco.Fill(reco_ene,mc_weif);
       fRecoND->Fill(reco_ene, mc_weif);
       cout <<"Found a tau in the ND MC!!!!!!!!!!" << reco_ene << endl;
    }
    for(int j = 0; j < size_elec; ++j){
       bin->GetMCNuEToNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       shower *= enu_shift;
       reco_ene = shower + mumom;
       sel_reco.Fill(reco_ene,mc_weif);
       fRecoND->Fill(reco_ene, mc_weif);
       cout <<"Found a CC-LIKE nue in the ND MC!!!!!!!!!!" << reco_ene << endl;
    }
    //end LLH added

  }//end loop over energy bins


  //for(uint i = 0; i < fweind.size(); ++i){

//     mc_weif       = fweind[i];
//     shower        = feshnd[i];
//     shower        = shower*enu_shift;
//     mumom         = femund[i];
//     reco_ene      = shower+mumom;
//     annpid        = fpidpnd[i];
//     cc_nc         = fccnd[i];
//     flavour       = fflavnd[i];
//     initial_state = finisnd[i];
//     true_enu      = fenutnd[i];

//     if(cc_nc == NCType::kCC
//        && flavour == 14
//        && initial_state <= 2){

//       selcc_true->Fill(true_enu,mc_weif);
//       recotruemat1->Fill(reco_ene,true_enu,mc_weif);
//       selcc_reco->Fill(reco_ene,mc_weif);
//     }

//     sel_reco->Fill(reco_ene,mc_weif);

//   }//end loop over near events

  double pur = 1.;
  double eff = 1.;
  double all = 0.;
  double selected = 0.;
  for(int i = 0; i < fnbins; i++){
    selected = selcc_reco.GetBinContent(i+1);
    all = sel_reco.GetBinContent(i+1);
    if(all > 0){
      pur = selected/all;
      purityn1->SetBinContent(i+1,pur);
    }

    selected = selcc_true.GetBinContent(i+1);
    all = fallccmc->GetBinContent(i+1);
    if(all > 0){
      eff = selected/all;
      efficiencyn1->SetBinContent(i+1,eff);
    }
  }
  //Efficiency check
  //---------------WHAT DOES THIS DO?
  double check = 1.;
  for(int i = 0; i < fnbins; i++){
    check=efficiencyn1->GetBinContent(i+1);
    if(check == 0 && i < fnbins)
      efficiencyn1->SetBinContent(i+1,efficiencyn1->GetBinContent(i+2));
  }

  Convert2DMatrixToProbability(recotruemat1);

  return;
}

//----------------------------------------------------------------------
//Why is this separate from GetEff??
void NCExtrapolationNS::RecoTrueEffPur(TH2F *truerecomat1, //cc sel cc
                                       TH2F *truerecomat2, //nc sel nc
                                       TH2F *truerecomat1b, //nc sel cc
                                       TH2F *truerecomat2b, //cc sel nc
                                       double enu_shift,
                                       int beam)
{
   cout <<"In RecoTrueEffPur! " << endl;

  truerecomat1->Reset();
  truerecomat2->Reset();

  truerecomat1b->Reset();
  truerecomat2b->Reset();

  double  mc_weif=1.;
  double  shower = 1.;
  double  mumom = 1.;
  double  reco_ene = 0.;
  int     flavour = 0;
  double  true_enu = 0.;
//   double   showernc = 1.;
//   int     annpid_f1p1 = 1;
//   int     cc_nc = 0;
//   int     initial_state = 0;

//   bool trcc   = false;
//   bool trcca  = false;
//   bool slcc   = false;
//   bool nslcc  = false;

//   bool evnc =false;
//   bool slnc =false;
//   bool nslnc=false;

  //size of the MC arrays
  int sizeSignal =  0;
  int sizeBG = 0;
  int sizeElec =0;

  //loop over bins in the beam
  NCEnergyBin *bin = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
    sizeSignal = bin->GetMCSignalVectorSize();
    sizeBG = bin->GetMCBackgroundVectorSize();
    sizeElec = bin->GetMCNuEToNuEVectorSize();

    //true cc selected as cc
    for(int j = 0; j < sizeSignal; ++j){
      bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      //only use nu_mu, no nu_mubar <= original case, LLH is including numubar!
      if(abs(flavour) == 14){
        truerecomat1->Fill(true_enu, reco_ene, mc_weif);
        fTrueFD->Fill(true_enu, mc_weif);
        fRecoFD->Fill(reco_ene, mc_weif);
        }
    }//end loop over signal events

    //true nc selected as cc
    for(int j = 0; j < sizeBG; ++j){
      bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      truerecomat2b->Fill(true_enu, reco_ene, mc_weif);
    }//end loop over background events

    //LLH adding
//    int size_tau = bin->GetMCNuMuToNuTauVectorSize(); //taus must be oscillated before adding, no osc on for debugging, therefore no taus!
//     for(int j = 0; j < size_tau; ++j){
//        bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
//        shower *= enu_shift;
//        reco_ene = shower + mumom;
//       fRecoFD->Fill(reco_ene, mc_weif);
//       cout <<"Finding FD tau!!!!" << endl;
//    }

    //true cc elec selected as cc
    for(int j = 0; j < sizeElec; ++j){
       bin->GetMCNuEToNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       shower *= enu_shift;
       reco_ene = shower + mumom;
       truerecomat2b->Fill(true_enu, reco_ene, mc_weif);
//        fRecoFD->Fill(reco_ene, mc_weif);
//        cout <<"Finding FD nue!!!!" << endl;
    }
    //end LLH added

  }//end loop over cc energy bins

  numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
    sizeSignal = bin->GetMCSignalVectorSize();
    sizeBG = bin->GetMCBackgroundVectorSize();
    sizeElec = bin->GetMCNuEToNuEVectorSize();

    //true nc selected as nc
    for(int j = 0; j < sizeSignal; ++j){
      bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      truerecomat2->Fill(true_enu,shower,mc_weif);
      fTrueFD_nc->Fill(true_enu, mc_weif);
      fRecoFD_nc->Fill(reco_ene, mc_weif);
    }//end loop over signal events

    //true cc selected as nc
    for(int j = 0; j < sizeBG; ++j){
      bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      truerecomat1b->Fill(true_enu,reco_ene,mc_weif);
    }//end loop over background events

    //true elec selected as nc, true cc only
    for(int j = 0; j < sizeElec; ++j){
       bin->GetMCNuEToNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       shower *= enu_shift;
       reco_ene = shower + mumom;
       truerecomat1b->Fill(true_enu, reco_ene, mc_weif);
    }

  }//end loop over nc energy bins

//   for(uint i = 0; i < fwei.size(); i++){
//     mc_weif       = fwei[i];
//     shower        = fesh[i];
//     shower        = shower*enu_shift;
//     showernc      = feshnc[i];
//     showernc      = showernc*enu_shift;
//     mumom         = femu[i];
//     reco_ene      = shower+mumom;
//     annpid_f1p1   = fpid_post[i];
//     cc_nc         = fccfd[i];
//     flavour       = fflavfd[i];
//     initial_state = finis[i];
//     true_enu      = fenut[i];

//     // CC Selection
//     trcc   = false;
//     trcca  = false;
//     slcc   = false;
//     nslcc  = false;

//     if(cc_nc == NCType::kCC
//        && initial_state <= 2
//        && flavour == 14)
//       trcc  = true;
//     if(cc_nc == NCType::kCC)
//       trcca = true;
//     if(annpid_f1p1 == NCType::kCC) slcc  = true;
//     if(annpid_f1p1 == NCType::kNC) nslcc = true;

//     // NC Selection
//     evnc =false;
//     slnc =false;
//     nslnc=false;

//     if(cc_nc == NCType::kNC) evnc = true;
//     if(annpid_f1p1 == NCType::kNC) slnc = true;
//     if(annpid_f1p1 == NCType::kCC) nslnc= true;

//     // CC selected as CC
//     if(trcc == true && slcc==true)
//       truerecomat1->Fill(true_enu,reco_ene,mc_weif);

//     // NC Selected as NC
//     if(evnc == true && slnc==true)
//       truerecomat2->Fill(true_enu,showernc,mc_weif);

//     //CC selected as NC
//     if(trcca==true && nslcc==true)
//       truerecomat1b->Fill(true_enu,showernc,mc_weif);

//     //NC selected as CC
//     if(evnc==true && nslnc==true)
//       truerecomat2b->Fill(true_enu,reco_ene,mc_weif);

//   }

  Convert2DMatrixToProbability(truerecomat1);
  Convert2DMatrixToProbability(truerecomat1b);
  Convert2DMatrixToProbability(truerecomat2);
  Convert2DMatrixToProbability(truerecomat2b);

  return;
}
//**********************************************************************
// TAUS RECO TO TRUE FOR BOTH CC AND ND NAMELY 2 HISTOS
void NCExtrapolationNS::RecoTrueEffPurTau(TH2F *truerecomat1,
                                          TH2F *truerecomat1b,
                                          double enu_shift,
                                          int beam)
{

  truerecomat1->Reset();
  truerecomat1b->Reset();

  double  mc_weif=1.;
  double   shower = 1.;
  double   mumom = 1.;
  double   reco_ene = 0.;
  int     flavour = 0;
  double   true_enu = 0.;
//   double   showernc = 1.;
//   int     annpid_f1p1 = 1;
//   int     cc_nc = 0;
//   int     initial_state = 0;
//   bool    selfl = false;
//   bool    selflnc = false;

  //size of the MC arrays
  int sizeTau =  0;

  //loop over bins in the beam
  //only true cc taus are included in the tau backgrounds as the nc are just nc events
  NCEnergyBin *bin = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();

    //true cc tau selected as cc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      reco_ene = shower + mumom;
      truerecomat1->Fill(true_enu,reco_ene,mc_weif);
    }//end loop over tau events
  }//end loop over cc energy bins

  numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();

    //true cc tau selected as nc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      shower *= enu_shift;
      truerecomat1b->Fill(true_enu,shower,mc_weif);
    }//end loop over tau events
  }//end loop over nc energy bins

//   for(uint i = 0; i < fweit.size(); i++){

//     mc_weif       = fweit[i];
//     shower        = fesht[i];
//     shower        = shower*enu_shift;
//     showernc      = feshnct[i];
//     showernc      = showernc*enu_shift;
//     mumom         = femut[i];
//     reco_ene      = shower+mumom;
//     annpid_f1p1   = fpidpt[i];
//     cc_nc         = fccfdt[i];
//     flavour       = fflavfdt[i];
//     initial_state = finist[i];
//     true_enu      = fenutt[i];

//     selfl = false;
//     selflnc = false;

//     // CC Selection
//     if(annpid_f1p1 == NCType::kCC) selfl =true;
//     // NC Selection
//     if(annpid_f1p1 == NCType::kNC) selflnc =true;
//     // CCTAUS selected as CC
//     if(cc_nc == NCType::kCC && flavour == 16){
//       if(selfl) truerecomat1->Fill(true_enu,reco_ene,mc_weif);
//     }
//     //CCTAUS selected as NC
//     if(cc_nc == NCType::kCC && flavour == 16){
//       if(selflnc)  truerecomat1b->Fill(true_enu,showernc,mc_weif);
//     }
//   }

  Convert2DMatrixToProbability(truerecomat1);
  Convert2DMatrixToProbability(truerecomat1b);

  return;
}

//********************************************************************************
void NCExtrapolationNS::GetEff(TH1F *effic_noosctrue, //cclike, cc (sig)
                               TH1F *effic_noosctrue_ncfit, //nclike, nc (sig)
                               TH1F *neffic_noosctrue, //nclike, cc (bg)
                               TH1F *neffic_noosctruenc, //cclike, nc (bg)
                               TH1F *ratccnc,
                               TH1F *rattau,
                               int beam)
{
  ratccnc->Reset();
  rattau->Reset();

  effic_noosctrue->Reset();
  effic_noosctrue_ncfit->Reset();

  neffic_noosctrue->Reset();
  neffic_noosctruenc->Reset();

//LLH deprecating  TH1F all("all",   "",fnbins,fbinvect);
  TH1F alltau("alltau","",fnbins,fbinvect);
//LLH deprecating  TH1F allnc("allnc", "",fnbins,fbinvect);

  double mc_weif = 1.;
  double shower = 0.;
  double mumom = 0.;
  int    flavour = 0;
  double true_enu = 0.;
//   double reco_ene = 0.;
//   int   annpid_f1p1 = 0;
//   int cc_nc = 0;
//   int initial_state = 0;

//   bool trcc =false;
//   bool trcca=false;
//   bool slcc =false;
//   bool nslcc=false;
//   bool evnc =false;
//   bool slnc =false;
//   bool nslnc=false;

  //size of the MC arrays
  int sizeTau =  0;
  int sizeSignal =  0;
  int sizeBG =  0;
  int sizeElec = 0;

  //loop over bins in the beam
  //only true cc taus are included in the tau backgrounds as the nc are just nc events
  NCEnergyBin *bin = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);  //Looping over CC selected
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
    sizeSignal = bin->GetMCSignalVectorSize();
    sizeBG = bin->GetMCBackgroundVectorSize();
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();
    sizeElec = bin->GetMCNuEToNuEVectorSize();

    //true cc selected as cc
    for(int j = 0; j < sizeSignal; ++j){
      bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      if(abs(flavour) == 14) effic_noosctrue->Fill(true_enu, mc_weif); //originally only numu, LLH including numubar!
    }//end loop over cc events

    //true nc selected as cc
    for(int j = 0; j < sizeBG; ++j){
      bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      neffic_noosctruenc->Fill(true_enu, mc_weif); //all 3 flavours!! (NC true only)
    }//end loop over nc events

    //true cc tau selected as cc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      alltau.Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
    }//end loop over tau events

    //true cc elec selected as cc
    for(int j = 0; j < sizeElec; ++j){
       bin->GetMCNuEToNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       neffic_noosctruenc->Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
    }//end loop over electron events
  }//end loop over cc energy bins

  numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);  //Looping over NC selected
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
    sizeSignal = bin->GetMCSignalVectorSize();
    sizeBG = bin->GetMCBackgroundVectorSize();
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();
    sizeElec = bin->GetMCNuEToNuEVectorSize();

    //true nc selected as nc
    for(int j = 0; j < sizeSignal; ++j){
      bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      effic_noosctrue_ncfit->Fill(true_enu, mc_weif); //numu and numubar (and neu neubar? - LLH)
    }//end loop over nc events

    //true cc selected as nc
    for(int j = 0; j < sizeBG; ++j){
      bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      neffic_noosctrue->Fill(true_enu, mc_weif);  //numu and numubar
    }//end loop over cc events

    //true tau selected as nc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      alltau.Fill(true_enu,mc_weif);  //fixme
    }//end loop over tau events

    //true elec selected as nc, true cc only
    for(int j = 0; j < sizeElec; ++j){
       bin->GetMCNuEToNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
       neffic_noosctrue->Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
    }//end loop over electron events
  }//end loop over nc energy bins

  rattau->Divide(&alltau,fallccmcf);
  ratccnc->Divide(fallncmcf,fallccmcf); //nc/cc - takes pred from allccmcf->allccncf, before eff corr

  feffic_numcheck = (TH1F*)effic_noosctrue->Clone("effic_numcheck"); //LLH adding, at this point should be identical to fTrueFD

  MSG("NCExtrapolationNS", Msg::kInfo) << " DONE WITH EFFICIENCIES " << endl;
  effic_noosctrue->Divide(fallccmcf);
  effic_noosctrue_ncfit->Divide(fallncmcf);

  neffic_noosctrue->Divide(fallccmcf);
  neffic_noosctruenc->Divide(fallncmcf);

}

//********************************************************
void NCExtrapolationNS::GetEfft(TH1F *effic_noosctrue,
                                TH1F *neffic_noosctrue,
                                int beam)
{
  effic_noosctrue->Reset();
  neffic_noosctrue->Reset();

  TH1F all("all",  "",fnbins,fbinvect);
  //TH1F *allnc   = new TH1F("allnc","",fnbins,fbinvect);

  double  mc_weif=1.;
  double  shower = 0.;
  double  mumom = 0.;
  int     flavour = 0;
  double  true_enu = 0.;
//   double  reco_ene = 0.;
//   int    annpid_f1p1 = 0;
//   int    cc_nc = 0;
//   int    initial_state = 0;

//   bool trcc = false;
//   bool slcc = false;
//   bool slnc = false;

  //size of the MC arrays
  int sizeTau =  0;

  //loop over bins in the beam
  //only true cc taus are included in the tau backgrounds as the nc are just nc events
  NCEnergyBin *bin = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();

    //true cc tau selected as cc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      all.Fill(true_enu,mc_weif);
      effic_noosctrue->Fill(true_enu,mc_weif);
    }//end loop over tau events
  }//end loop over cc energy bins

  numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
    sizeTau = bin->GetMCNuMuToNuTauVectorSize();

    //true cc tau selected as cc
    for(int j = 0; j < sizeTau; ++j){
      bin->GetMCNuMuToNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
      all.Fill(true_enu,mc_weif);
      neffic_noosctrue->Fill(true_enu,mc_weif);
    }//end loop over tau events
  }//end loop over cc energy bins

//   for(uint i = 0; i < fweit.size(); i++){

//     mc_weif       = fweit[i];
//     mc_weif       = fweit[i];
//     shower        = fesht[i];
//     mumom         = femut[i];
//     reco_ene      = shower+mumom;

//     annpid_f1p1   = fpidpt[i];
//     cc_nc         = fccfdt[i];
//     flavour       = fflavfdt[i];
//     initial_state = finist[i];
//     true_enu      = fenutt[i];

//     trcc = false;
//     slcc = false;
//     slnc = false;

//     if(cc_nc == NCType::kCC && flavour == 16) trcc = true;
//     if(annpid_f1p1 == NCType::kCC) slcc = true;
//     if(annpid_f1p1 == NCType::kNC) slnc = true;

//     if(trcc==true)                all->Fill(true_enu,mc_weif);
//     if(slcc==true  && trcc==true) effic_noosctrue->Fill(true_enu,mc_weif);
//     if(slnc==true  && trcc==true) neffic_noosctrue->Fill(true_enu,mc_weif);
//   }

  MSG("NCExtrapolationNS", Msg::kInfo) << " DONE WITH EFFICIENCIES TAU " << endl;
  effic_noosctrue->Divide(&all);
  neffic_noosctrue->Divide(&all);

  return;
}

//***************************************************************

//LLH Note: this only works if inihist has smaller bins that fit exactly inside of rebinhist bins
void NCExtrapolationNS::Rebindata(TH1F *inihist,
                                  TH1F &rebinhist)
{
  //loop over the output histogram and find the edges of the
  //bins.  loop over the input histogram and put all bins that land
  //within those edges into the output histogram
  double binCenter = 0.;
  double binContent = 0.;
  double outBinLow = 0.;
  double outBinHigh = 0.;
  double inBinLow = 0.;
  double inBinHigh = 0.;
  rebinhist.Reset();
  for(int i = 0; i < rebinhist.GetXaxis()->GetNbins(); ++i){
    outBinLow = rebinhist.GetXaxis()->GetBinLowEdge(i+1);
    outBinHigh = rebinhist.GetXaxis()->GetBinUpEdge(i+1);
    for(int j = 0; j < inihist->GetXaxis()->GetNbins(); ++j){
      binContent = inihist->GetBinContent(j+1);
      binCenter = inihist->GetBinCenter(j+1);
      inBinLow = inihist->GetXaxis()->GetBinLowEdge(j+1);
      inBinHigh = inihist->GetXaxis()->GetBinUpEdge(j+1);

      //if the input histogram bin is withint the output bin,
      //add the content to the output histogram
      if(binCenter >= outBinLow && binCenter <= outBinHigh){
        if(inBinLow >= outBinLow && inBinHigh <= outBinHigh){
          rebinhist.Fill(binCenter, binContent);
          MSG("NCExtrapolationNS", Msg::kDebug)
            << rebinhist.GetName() << " " << outBinLow
            << " " << outBinHigh << " " << inBinLow
            << " " << inBinHigh << " " << rebinhist.GetBinContent(i+1)
            << " " << binContent << endl;
        }
        else if(inBinLow >= outBinLow && inBinHigh >= outBinHigh)
          MSG("NCExtrapolationNS", Msg::kWarning) << "input bin goes over output bin high edge "
                                                  << binCenter << " "
                                                  << outBinLow << "/" << inBinLow << " "
                                                  << inBinHigh << "/" << outBinHigh << endl;
        else if(inBinLow <= outBinLow && inBinHigh <= outBinHigh)
          MSG("NCExtrapolationNS", Msg::kWarning) << "input bin goes under output bin low edge "
                                                  << binCenter << " "
                                                  << outBinLow << "/" << inBinLow << " "
                                                  << inBinHigh << "/" << outBinHigh << endl;
      }//end if input bin center is within output bin

    }//end loop over input histogram
  }//end loop over output histogram

  return;
}

//----------------------------------------------------------------------
//Rebindata now works for any input histogram as it uses the information
//from the rebinhist and inihist to determine how to combine the bins in
//inihist
void NCExtrapolationNS::Rebin(TH1F *inihist,
                              TH1F &rebinhist)
{
  Rebindata(inihist, rebinhist);

//   int kkk=1;

//   for(int kk=1;kk<200;kk++){
//     if(kk==3){
//       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1)+inihist ->GetBinContent(kk-2);
//       rebinhist.SetBinContent(kkk,reco2);
//     }
//     if(kk>3 && kk<=39 && kk%2==1){
//       kkk=kkk+1;
//       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1);
//       rebinhist.SetBinContent(kkk,reco2);
//     }
//     if(kk>39 && kk%4==3 && kk<=119){ // 75
//       kkk=kkk+1;
//       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1)+inihist->GetBinContent(kk-2)+inihist ->GetBinContent(kk-3);
//       rebinhist.SetBinContent(kkk,reco2);
//     }

//   }

  return;
}

//----------------------------------------------------------------------
float NCExtrapolationNS::chisq(TH1F *obs,
                               TH1F *exp,
                               TH1F *expnc)
{

  // calculate chisq between 2 histograms - use Poisson form

  float chisquare=0;
  float nobs = 0.;
  float nexp = 0.;
  float nexpnc = 0.;

  for(int i = 0; i < obs->GetNbinsX(); i++){
    nobs = obs->GetBinContent(i+1);
    nexp = exp->GetBinContent(i+1);

    nexpnc = expnc->GetBinContent(i+1);

    if(nexp > 0){
      // LOG LIKELIHOOD
      if(nobs > 0)  chisquare += 2*(nexp-nobs)+2*nobs*log(nobs/nexp);
      if(nobs == 0) chisquare += 2*(nexp-nobs);
    }

// Temp - LLH
//     MSG("NCExtrapolationNS", Msg::kInfo) << "find chi^2 " << obs->GetBinCenter(i+1)
//                                       << " " << nobs << " / " << nexp
//                                       << " " << nexp/nobs << " "
//                                       << " " << chisquare << endl;
  }

  return chisquare;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::truetoreco(TH1F *truth,
                                   TH1F& reco,
                                   TH2F *matrix)
{
  //Adding this back (got lost in translation for some reason) - LLH
   reco.Reset();

  // use matrix to convert from true to reco
  int nbins = truth->GetNbinsX();

  float count = 0.;
  for(int i = 0; i < nbins; i++){
    count = 0;
    for(int j = 0; j < nbins; j++){
      count += truth->GetBinContent(j+1)*matrix->GetBinContent(j+1, i+1);
    }
    reco.SetBinContent(i+1,count);

// Temp
//     MAXMSG("NCExtrapolationNS", Msg::kInfo, nbins) << "far true to reco "
//                                                 << ":  bin_center = " << reco.GetBinCenter(i+1) << " "
//                                                 << ", reco = " << count << ", truth = " << truth->GetBinContent(i+1)
//                                                 << endl;
  }

  return;
}

//----------------------------------------------------------------------
//this method uses the value of sin^2(1.27Delta m^2 L/E) appropriate for each bin
//in the provided vector.
void NCExtrapolationNS::oscweight(std::vector<double> &sinSqrDeltaMSqr,
                                  double sinSqr2Theta, double fSterile,
                                  TH1F *noosc, TH1F *osc, int flavor,
                                  int interactionType)
{
  double prob = 1.;
  double noOscVal = 0.;

  for(uint i = 0; i < sinSqrDeltaMSqr.size(); ++i){
    noOscVal = noosc->GetBinContent(i+1);

    if(TMath::Abs(flavor) == 14){
      if(interactionType == NCType::kCC)
        prob = 1. - sinSqr2Theta*sinSqrDeltaMSqr[i];
      else if(interactionType == NCType::kNC)
        prob = 1. - sinSqr2Theta*fSterile*sinSqrDeltaMSqr[i];
    }
    else if(TMath::Abs(flavor) == 16){
      if(interactionType == NCType::kCC)
        prob = (1. - fSterile)*sinSqr2Theta*sinSqrDeltaMSqr[i];
      else if(interactionType == NCType::kNC)
        prob = 0.;
    }

    osc->SetBinContent(i+1, prob*noOscVal);
  }//end loop over energy bins

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::oscweightOld(TH1F *noosc,
                                     TH1F &osc,
                                     float dm,
                                     float st,
                                     int method,
                                     int what)
{

  // what = 1 cc disappearance
  // what = 2 nc disappearance
  // what = 11 taus appearance
  // what = -1 decay

  // reweight spectrum with oscillation parameters
  // interpolation method:  0 - none,  1 - quadratic,  anything else - linear
  //osc.Reset();
  float distance = NCType::kBaseLineFar;

  int nbins = noosc->GetNbinsX();
  float prob = 0.;
  float binsize = 0.25;
  float wgt1 = 0.;
  float wgt2 = 0.;
  float wgt3 = 0.;
  float de = 0.;
  float de2 = 0.;

  float avgosc = 0.;
  float avgosc2 = 0.;
  float avgwgt1 = 0.;
  float avgwgt2 = 0.;
  float enu = 0.;
  float delta = 0.;
  float linwgt = 0.;
  float quadwgt = 0.;

  for (int i = 0; i < nbins; i++) {

    if(what == 1 || what == 0)
      prob = 1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/noosc->GetBinCenter(i+1)),2);
    if(what == 11)
      prob = st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/noosc->GetBinCenter(i+1)),2);
    if(what == -1)
      prob=TMath::Power(st+(1-st)*TMath::Exp(-dm*distance/(2*noosc->GetBinCenter(i+1))),2);
    if(what == -11)
      prob = 0.;

    binsize = 0.25;
    if (i == 0)    binsize = 0.5;
    if (i == 199)  binsize = 200.-49.75;

    if(i > nbins-2 || method == 0){
      osc.SetBinContent(i+1, noosc->GetBinContent(i+1)*prob);
    }
    else{

      wgt1 = noosc->GetBinContent(i+1);
      wgt2 = noosc->GetBinContent(i+2);
      wgt3 = noosc->GetBinContent(i+3);

      if(i == 0){
        wgt2 = noosc->GetBinContent(i+2)+noosc->GetBinContent(i+3);
        wgt3 = noosc->GetBinContent(i+4)+noosc->GetBinContent(i+5);

      }

      de = wgt2-wgt1;
      de2 = wgt3-2*wgt2+wgt1;

      avgosc = 0.;
      avgosc2 = 0.;
      avgwgt1 = 0.;
      avgwgt2 = 0.;

      for(int k = 0; k < 10; k++){

        enu = noosc->GetBinLowEdge(i+1)+k*binsize/10.+binsize/20.;

        delta = enu-noosc->GetBinCenter(i+1);

        linwgt = wgt1+de*delta;
        quadwgt = wgt1+de*delta+de2/2*TMath::Power(delta,2);

        if(linwgt < 0) linwgt = 0.;
        if(quadwgt < 0) quadwgt = 0.;

        if(what == 1 || what == 0){
          avgosc += quadwgt*(1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
          avgosc2 += linwgt*(1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
        }
        if(what==11){
          avgosc += quadwgt*(st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
          avgosc2 += linwgt*(st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
        }
        if(what==-1){
          avgosc += quadwgt*TMath::Power(st+(1-st)*exp(-dm*distance/(2*enu)),2);
          avgosc2 += linwgt*TMath::Power(st+(1-st)*exp(-dm*distance/(2*enu)),2);
        }
        if(what==-11){
          avgosc += quadwgt*0.;
          avgosc2 += linwgt*0.;
        }
        avgwgt1 += quadwgt;
        avgwgt2 += linwgt;
      }

      if(method == 1 && avgosc > 0){
        if(avgwgt1 > 0) avgosc /= avgwgt1;
        else{
          MSG("NCExtrapolationNS", Msg::kWarning) << " attempted divide by 0 to get average "
                                                  << " oscillation for bin -  not averaged "
                                                  << avgosc << " " << avgwgt1 << " "
                                                  << osc.GetBinCenter(i+1) << " "
                                                  << noosc->GetBinContent(i+1) << endl;
        }
        osc.SetBinContent(i+1, avgosc*noosc->GetBinContent(i+1));
      }
      else if(avgwgt2 > 0){
        if(avgwgt2 > 0) avgosc2 /= avgwgt2;
        else{
          MSG("NCExtrapolationNS", Msg::kWarning) << " attempted divide by 0 to get average "
                                                  << " oscillation for bin -  not averaged "
                                                  << avgosc2 << " " << avgwgt2 << endl;
        }
        osc.SetBinContent(i+1, avgosc2*noosc->GetBinContent(i+1));
      }

    }
  }//end loop over energy bins

  if((st == 0 || dm == 0) && (what == 1 || what == 0)){
    for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,noosc->GetBinContent(i+1));
  }

  if((st == 0 || dm == 0) && (what == 11 || what == -11)){
    for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,0.);
  }

  if(dm == 0 && what == -1){
    for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,noosc->GetBinContent(i+1));
  }

  return;
}

//----------------------------------------------------------------------
//matricies are either reco to true or true to reco.  if reco to true,
//the x axis is reconstructed energy, the y axis is true energy.
//if true to reco, the x axis is true energy, the y axis is reco'ed
//energy
void NCExtrapolationNS::Convert2DMatrixToProbability(TH2F *matrix)
{
  double sum[fnbins];
  double maxf[fnbins][fnbins];

  //initialize
  for(int j = 0; j < fnbins; j++){
    sum[j] = 0.;
    for(int i = 0; i < fnbins; i++) maxf[i][j] = 0.;
  }

  //Make True to reco or reco to true
  //this step takes the y bins and puts them in the
  //x bins and vice versa
  for(int i = 0; i < fnbins; i++){
    for(int j = 0; j < fnbins; j++){
      maxf[j][i] = matrix->GetBinContent(i+1, j+1);
    }
  }

  // MAKE 2D MATRIX A PROBABILITY ONE
//LLH in here twice  for(int i = 0; i < fnbins; i++) sum[i] = 0.;

  //sum over the new x bins for each y value
  for(int j = 0; j < fnbins; j++){
    for(int i = 0; i < fnbins; i++){
      sum[j] += maxf[i][j];
    }
  }

  //each x,y bin is the fraction of total for
  //given row in y.  set bin content in matrix to
  //be true,reco (reco,true)
  for(int i = 0; i < fnbins; i++){
    for(int j = 0; j < fnbins; j++){
      if(sum[j] > 0.){
        maxf[i][j] /= sum[j];
        matrix->SetBinContent(j+1, i+1, maxf[i][j]);
      }
    }
  }

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationNS::FillDataHistogramFromEnergyBins(int beam,
                                                        int nccc,
                                                        TH1F *hist)
{
  hist->Reset();

  NCEnergyBin *bin = 0;
  int sizeData = 0;
  int numBins = fBeams[beam].GetNumberEnergyBins(nccc);
  double energy = 0.;
  double weight = 0.;
  for(int i = 0; i < numBins; ++i){
    bin = fBeams[beam].GetEnergyBin(i, nccc);
    sizeData = bin->GetDataVectorSize();
    for(int j = 0; j < sizeData; ++j){
      bin->GetDataInformation(energy, weight, j);
      hist->Fill(energy, weight);
    }
  }//end loop over bins

  return;
}

---