NCExtrapolationBeamMatrix.cxx

Go to the documentation of this file.

//
// NCExtrapolationBeamMatrix.cxx
//
//
// R. Pittam 10/2008
//
// $Id: NCExtrapolationBeamMatrix.cxx,v 1.7 2009/12/10 17:07:06 gmieg Exp $

#include "NCUtils/Extrapolation/NCExtrapolationBeamMatrix.h"
#include "NCUtils/Extrapolation/NCEnergyBin.h"
#include "NCUtils/Extrapolation/NCPOTCounter.h"
#include "NCUtils/NCUtility.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 "TCanvas.h"
#include "TDirectory.h"
#include "TGraphErrors.h"
#include "TMath.h"
#include "TMinuit.h"

#include <fstream>
#include <cassert>
#include <cmath>

CVSID("$Id: NCExtrapolationBeamMatrix.cxx,v 1.7 2009/12/10 17:07:06 gmieg Exp $");

#include "NCExtrapolationModule.h"
#include "NCSpectrumInterpolator.h"

REGISTER_NCEXTRAPOLATION(NCExtrapolationBeamMatrix, BeamMatrix)

NCExtrapolationBeamMatrix* NCExtrapolationBeamMatrix::bmFit = NULL;

using namespace NC::Utility;

const int kNumTrueEnergyBins = 1000;

enum EBeamMatrixEventType{
  kNCSig     = 0, //NC spectrum
  kCCSig     = 1, //CC spectrum
  kNCBkg     = 2, //NC background in NC spectrum
  kCCBkg     = 3, //CC background in CC spectrum
  kNCTau     = 4, //NC tau background in NC spectrum
  kCCTau     = 5, //CC tau background in CC spectrum
  kNCBeamNue = 6, //NC nue background in NC spectrum
  kCCBeamNue = 7, //CC nue background in CC spectrum
  kNCOscNue  = 8, //NC oscillated nue background in NC spectrum
  kCCOscNue  = 9, //CC oscillated nue background in CC spectrum
  kNCWrongSignBkg  = 10, //NC wrong sign background in NC spectrum
  kCCWrongSignBkg  = 11 //CC wrong sign background in CC spectrum
};

enum EBeamMatrixHistoType{
  kNCRecoEnergy_ND  = 0, // NC Histogram for Reco Spectrum in the ND
  kCCRecoEnergy_ND  = 1, // CC Histogram for Reco Spectrum in the ND
  kNCRecoBkg_ND  = 2, // NC Background in Reco Energy in the ND
  kCCRecoBkg_ND  = 3, // CC Background in Reco Energy in the ND
  kNCPurity_ND  = 4, // NC Histogram for Purity correction in the ND
  kCCPurity_ND  = 5, // CC Histogram for Purity correction in the ND
  kNCRecoEnergyPurCorr_ND  = 6, // NC Histogram for Purity corrected Reco Spectrum in the ND
  kCCRecoEnergyPurCorr_ND  = 7, // CC Histogram for Purity corrected Reco Spectrum in the ND
  kNCTrueEnergyPurCorr_ND  = 8, // NC Histogram for Purity corrected True Spectrum in the ND
  kCCTrueEnergyPurCorr_ND  = 9, // CC Histogram for Purity corrected True Spectrum in the ND
  kNCTrueSelEfficiency_ND  = 10, // NC Histogram for Selection Efficency correction in the ND (done in true energy)
  kCCTrueSelEfficiency_ND  = 11, // CC Histogram for Selection Efficency correction in the ND (done in true energy)
  kNCTrueEnergySelEffCorr_ND  = 12, // NC Histogram for Purity and Selection Efficiency corrected True Spectrum in the ND
  kCCTrueEnergySelEffCorr_ND  = 13, // CC Histogram for Purity and Selection Efficiency corrected True Spectrum in the ND
  kNCTrueRecoEfficiency_ND  = 14, // NC Histogram for Reconstruction Efficiency correction in the ND (done in true energy)
  kCCTrueRecoEfficiency_ND  = 15, // CC Histogram for Reconstruction Efficency correction in the ND (done in true energy)
  kNCTrueEnergyEffCorr_ND  = 16, // NC Histogram for Purity, Selection & Reconstruction Efficiency correction in the ND (done in true energy)
  kCCTrueEnergyEffCorr_ND  = 17, // CC Histogram for Purity, Selection & Reconstruction Efficiency correction in the ND (done in true energy)
  kNCTrueEnergy_ND  = 18, // NC Histogram for the ACTUAL (From MC) Purity and Efficiency corrected True Spectrum in the ND needed to make efficiency histogram
  kCCTrueEnergy_ND  = 19, // CC Histogram for the ACTUAL (From MC) Purity and Efficiency corrected True Spectrum in the ND needed to make efficiency histogram
  kNCTrueEnergyBkg_ND  = 20, // NC Histogram for the Background from MC   
  kCCTrueEnergyBkg_ND  = 21, // CC Histogram for the Background from MC
  kNCTrueEnergyFlux_ND  = 22, // NC Histogram for the ND Neutrino Flux after pot mass and x-section corrections
  kCCTrueEnergyFlux_ND  = 23, // CC Histogram for the ND Neutrino Flux after pot mass and x-section corrections
  kNCTrueEnergyMatrix_FD  = 24, // NC Histogram for the FD prediction after the Beam Matrix
  kCCTrueEnergyMatrix_FD  = 25, // CC Histogram for the FD prediction aftet the Beam Matrix
  kNCTrueEnergyFlux_FD  = 26, // NC Histogram for the FD Neutrino Flux after pot mass and x-section corrections
  kCCTrueEnergyFlux_FD  = 27, // CC Histogram for the FD Neutrino Flux after pot mass and x-section corrections
  kNCTrueSelEfficiency_FD  = 28, // NC Histogram for Selection Efficency correction in the ND (done in true energy)
  kCCTrueSelEfficiency_FD  = 29, // CC Histogram for Selection Efficency correction in the ND (done in true energy)
  kNCTrueEnergySelEffCorr_FD  = 30, // NC Histogram for Purity and Selection Efficiency corrected True Spectrum in the ND
  kCCTrueEnergySelEffCorr_FD  = 31, // CC Histogram for Purity and Selection Efficiency corrected True Spectrum in the ND
  kNCTrueRecoEfficiency_FD  = 32, // NC Histogram for Reconstruction Efficiency correction in the ND (done in true energy)
  kCCTrueRecoEfficiency_FD  = 33, // CC Histogram for Reconstruction Efficiency correction in the ND (done in true energy)
  kNCTrueEnergyEffCorr_FD  = 34, // NC Histogram for Purity and Efficiency corrected True Spectrum in the ND
  kCCTrueEnergyEffCorr_FD  = 35, // CC Histogram for Purity and Efficiency corrected True Spectrum in the ND
  kNCTrueEnergy_FD  = 36, // NC Histogram for the ACTUAL Purity and Efficiency corrected True Spectrum in the FD used as nominal MC
  kCCTrueEnergy_FD  = 37, // CC Histogram for the ACTUAL Purity and Efficiency corrected True Spectrum in the FD used as nominal MC
  kNCTrueEnergyXsecFit_FD  = 38, // NC Histogram for the ACTUAL Purity and Efficiency corrected True Spectrum in the FD needed normalise FD TrueToReco
  kCCTrueEnergyXsecFit_FD  = 39, // CC Histogram for the ACTUAL Purity and Efficiency corrected True Spectrum in the FD needed normalise FD TrueToReco
  kNCTrueEnergyBkg_FD  = 40, // NC Histogram for the Background to make above 2 histograms needed to make efficiency histogram 
  kCCTrueEnergyBkg_FD  = 41, // CC Histogram for the Background to make above 2 histograms needed to make efficiency histogram
  kNCRecoEnergy_FD = 42, // NC Histogram for Reco Spectrum in the FD
  kCCRecoEnergy_FD = 43, // CC Histogram for Reco Spectrum in the FD
  kNCRecoEnergyPurCorr_FD  = 44, // NC Histogram for Purity corrected Reco Spectrum in the FD
  kCCRecoEnergyPurCorr_FD  = 45, // CC Histogram for Purity corrected Reco Spectrum in the FD
  kNCPurity_FD  = 46, // NC Histogram for Purity correction in the FD
  kCCPurity_FD  = 47, // CC Histogram for Purity correction in the FD
  kNCRecoBkg_FD  = 48, // NC Histogram for background in reco energy
  kCCRecoBkg_FD  = 49, // CC Histogram for background in reco energy
  kNCRecoEnergyPred_FD = 50, // NC BM prediction in the FD
  kCCRecoEnergyPred_FD = 51, // CC BM prediction in the FD
  kNCRecoEnergyBeamNueBkg_FD = 52, // CC Histogram for Reco Spectrum in the FD
  kCCRecoEnergyBeamNueBkg_FD = 53, // CC Histogram for Reco Spectrum in the FD
  kNCRecoEnergyBeamNueBkg_ND = 54, // CC Histogram for Reco Spectrum in the ND
  kCCRecoEnergyBeamNueBkg_ND = 55, // CC Histogram for Reco Spectrum in the ND
  kNCRecoEnergyTau_ND = 56, // CC Histogram for Reco Spectrum in the ND
  kCCRecoEnergyTau_ND = 57, // CC Histogram for Reco Spectrum in the ND
  kNCRecoEnergyOscNue_ND = 58, // CC Histogram for Reco Spectrum in the ND
  kCCRecoEnergyOscNue_ND = 59, // CC Histogram for Reco Spectrum in the ND
  kNCSelRecoEnergy_ND  = 60, // NC Histogram for Reco Spectrum in the ND
  kCCSelRecoEnergy_ND  = 61, // CC Histogram for Reco Spectrum in the ND
  kNCSelRecoEnergy_FD  = 62, // NC Histogram for Reco Spectrum in the FD
  kCCSelRecoEnergy_FD  = 63 // CC Histogram for Reco Spectrum in the FD
};

enum EBeamMatrixDataHistoType{
  kNCDataRecoEnergy_ND  = 0, //NC Histogram for Reco Spectrum in the ND
  kCCDataRecoEnergy_ND  = 1, //CC Histogram for Reco Spectrum in the ND
  kNCDataRecoEnergyPurCorr_ND  = 2, //NC Histogram for Purity corrected Reco Spectrum in the ND
  kCCDataRecoEnergyPurCorr_ND  = 3, //CC Histogram for Purity corrected Reco Spectrum in the ND
  kNCDataTrueEnergyPurCorr_ND  = 4, //NC Histogram for Purity corrected True Spectrum in the ND
  kCCDataTrueEnergyPurCorr_ND  = 5, //CC Histogram for Purity corrected True Spectrum in the ND
  kNCDataTrueEnergySelEffCorr_ND  = 6, //NC Histogram for Purity and  selection Efficiency corrected True Spectrum in the ND
  kCCDataTrueEnergySelEffCorr_ND  = 7, //CC Histogram for Purity and selection Efficiency corrected True Spectrum in the ND
  kNCDataTrueEnergyEffCorr_ND  = 8, //NC Histogram for Purity and selection Efficiency corrected True Spectrum in the ND
  kCCDataTrueEnergyEffCorr_ND  = 9, //CC Histogram for Purity and selection Efficiency corrected True Spectrum in the ND
  kNCDataTrueEnergyFlux_ND  = 10, // NC Histogram for the ND Flux after pot mass and x-section corrections
  kCCDataTrueEnergyFlux_ND  = 11, // CC Histogram for the ND Flux after pot mass and x-section corrections
  kNCDataTrueEnergyMatrix_FD  = 12, // NC Histogram for the FD prediction after the Beam Matrix
  kCCDataTrueEnergyMatrix_FD  = 13, // CC Histogram for the FD prediction aftet the Beam Matrix
  kNCDataTrueEnergyFlux_FD  = 14, // NC Histogram for the FD Flux after pot mass and x-section corrections
  kCCDataTrueEnergyFlux_FD  = 15, // CC Histogram for the FD Flux aftet pot mass and x-section corrections
  kNCDataTrueEnergySelEffCorr_FD  = 16, //NC Histogram for Purity and selection Efficiency corrected True Spectrum in the ND
  kCCDataTrueEnergySelEffCorr_FD  = 17, //CC Histogram for Purity and selection Efficiency corrected True Spectrum in the ND
  kNCDataTrueEnergyEffCorr_FD  = 18, //NC Histogram for Purity and selection and reconstruction Efficiency corrected True Spectrum in the ND
  kCCDataTrueEnergyEffCorr_FD  = 19, //CC Histogram for Purity and selection and reconstruction Efficiency corrected True Spectrum in the ND
  kNCDataRecoEnergyPurCorr_FD = 20, // NC Histogram for Purity corrected Reco Spectrum in the FD
  kCCDataRecoEnergyPurCorr_FD = 21, //CC Histogram for Purity corrected Reco Spectrum in the FD
  kNCDataRecoEnergyPred_FD = 22,  // NC BM prediction in the FD
  kCCDataRecoEnergyPred_FD = 23, // CC BM prediction in the FD
  kNCDataRecoEnergySig_FD = 24, // NC BM predicted signal only part all the others can be got from the canvases
  kCCDataRecoEnergySig_FD = 25, // CC BM predicted signal only part all the others can be got from the canvases
  kNCDataRecoEnergyWrongSignBkg_FD = 26, // NC BM predicted wrong sign bkg all the others can be got from the canvases
  kCCDataRecoEnergyWrongSignBkg_FD = 27 // CC BM predicted wrong sign bkg all the others can be got from the canvases
};

enum EBeamMatrixTruthHistoType{
  kNCRecoTrueSig_ND = 0, // NC Reco to true matrix flavour corrected in the ND
  kCCRecoTrueSig_ND = 1, // CC Reco to true matrix flavour corrected in the ND
  kNCRecoTrueBkg_ND = 2, // NC Reco to true matrix flavour corrected wrong type in the ND
  kCCRecoTrueBkg_ND = 3, // CC Reco to true matrix flavour corrected wrong type in the ND
  kNCRecoTrueNorm_ND = 4, // NC Reco to true matrix Normalised in the ND
  kCCRecoTrueNorm_ND = 5, // CC Reco to true matrix Normalised in the ND

  kNCTrueRecoSig_FD = 6, // NC Reco to true matrix flavour corrected in the FD
  kCCTrueRecoSig_FD = 7, // CC Reco to true matrix flavour corrected in the FD
  kNCTrueRecoBkg_FD = 8, // NC Reco to true matrix flavour corrected wrong type in the FD
  kCCTrueRecoBkg_FD = 9, // CC Reco to true matrix flavour corrected wrong type in the FD
  kNCTrueRecoNorm_FD = 10, // NC Reco to true matrix flavour corrected in the FD
  kCCTrueRecoNorm_FD = 11 // CC Reco to true matrix flavour corrected in the FD
};

enum EBeamMatrixFitXsecType{
  kNCDataLE = 0, 
  kNCMCLE = 1, 
  kNCDataME = 2, 
  kNCMCME = 3, 
  kNCDataHE = 4, 
  kNCMCHE = 5, 
  kCCDataLE = 6, 
  kCCMCLE = 7, 
  kCCDataME = 8, 
  kCCMCME = 9, 
  kCCDataHE = 10, 
  kCCMCHE = 11 
};

int bmeqold = -1 ;

//......................................................................
NCExtrapolationBeamMatrix::NCExtrapolationBeamMatrix()
{
  fNCalls = 0;
  fDoExtrapolation = true;
  
  int numBins = kNumEnergyBinsFar;
  double bins[kNumEnergyBinsFar+1];
  
  
  //use NCType to fill array of bins as it has the right binning already
  for( int i = 0; i < numBins+1; ++i ){
    bins[i] = kEnergyBinsFar[i];
  }
  
  // Create a true energy binning scheme for the
  // FD prediction to oscillate
  
  //1000
  for( int i = 0; i < 200; ++i )
    true_bins.push_back(i*0.01);
  for( int i = 80; i < 801; ++i )
    true_bins.push_back(i*0.025);
  for( int i = 41; i < 81; ++i )
    true_bins.push_back(i*0.5);
  for( int i = 42; i < 120; i+=2 )
    true_bins.push_back(i);
  true_bins.push_back(120) ;
  
  
  // needed for ND fitting
  std::vector< TH1D* > NDByEnergyForXsecLE; 
  NDByEnergyForXsecLE.push_back(new TH1D(" NC0to4LE" , " NC0to4LE" , numBins, bins) );
  NDByEnergyForXsecLE.push_back(new TH1D(" NC5to8LE" , " NC5to8LE" , numBins, bins) );
  NDByEnergyForXsecLE.push_back(new TH1D(" NC8to15LE" , "NC8to15LE" , numBins, bins) );
  NDByEnergyForXsecLE.push_back(new TH1D(" NCover15LE" , "NCover15LE" , numBins, bins) );
  NDByEnergyForXsecLE.push_back(new TH1D(" NCnotMuLE" , "NCnotMuLE" , numBins, bins) );

  std::vector< TH1D* > NDByEnergyForXsecME; 
  NDByEnergyForXsecME.push_back(new TH1D(" NC0to4ME" , " NC0to4ME" , numBins, bins) );
  NDByEnergyForXsecME.push_back(new TH1D(" NC5to8ME" , " NC5to8ME" , numBins, bins) );
  NDByEnergyForXsecME.push_back(new TH1D(" NC8to15ME" , "NC8to15ME" , numBins, bins) );
  NDByEnergyForXsecME.push_back(new TH1D(" NCover15ME" , "NCover15ME" , numBins, bins) );
  NDByEnergyForXsecME.push_back(new TH1D(" NCnotMuME" , "NCnotMuME" , numBins, bins) );
  
  std::vector< TH1D* > NDByEnergyForXsecHE; 
  NDByEnergyForXsecHE.push_back(new TH1D(" NC0to4HE" , " NC0to4HE" , numBins, bins) );
  NDByEnergyForXsecHE.push_back(new TH1D(" NC5to8HE" , " NC5to8HE" , numBins, bins) );
  NDByEnergyForXsecHE.push_back(new TH1D(" NC8to15HE" , "NC8to15HE" , numBins, bins) );
  NDByEnergyForXsecHE.push_back(new TH1D(" NCover15HE" , "NCover15HE" , numBins, bins) );
  NDByEnergyForXsecHE.push_back(new TH1D(" NCnotMuHE" , "NCnotMuHE" , numBins, bins) );
  
  ND_XsecByParams.push_back(NDByEnergyForXsecLE);
  ND_XsecByParams.push_back(NDByEnergyForXsecME);
  ND_XsecByParams.push_back(NDByEnergyForXsecHE);


  ND_XSectionFitRew.push_back(new TH1D("NCDataLE" , "NCDataLE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("NCMCLE" , "NCMCLE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("NCDataME" , "NCDataME" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("NCMCME" , "NCMCME" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("NCDataHE" , "NCDataHE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("NCMCHE" , "NCMCHE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCDataLE" , "CCDataLE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCMCLE" , "CCMCLE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCDataME" , "CCDataME" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCMCME" , "CCMCME" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCDataHE" , "CCDataHE" , numBins, bins) );
  ND_XSectionFitRew.push_back(new TH1D("CCMCHE" , "CCMCHE" , numBins, bins) );

  ND_XSectionFit.push_back(new TH1D("NCDataLEOrig" , "NCDataLEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("NCMCLEOrig" , "NCMCLEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("NCDataMEOrig" , "NCDataMEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("NCMCMEOrig" , "NCMCMEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("NCDataHEOrig" , "NCDataHEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("NCMCHEOrig" , "NCMCHEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCDataLEOrig" , "CCDataLEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCMCLEOrig" , "CCMCLEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCDataMEOrig" , "CCDataMEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCMCMEOrig" , "CCMCMEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCDataHEOrig" , "CCDataHEOrig" , numBins, bins) );
  ND_XSectionFit.push_back(new TH1D("CCMCHEOrig" , "CCMCHEOrig" , numBins, bins) );


  dataNDhistLE = new TH1D("dataNDhistLE","dataNDhistLE", numBins, bins);
  dataNDhistME = new TH1D("dataNDhistME","dataNDhistME", numBins, bins);
  dataNDhistHE = new TH1D("dataNDhistHE","dataNDhistHE", numBins, bins);
  
  //Efficency histograms read from file
  selEffNDHist =  new TH1D();
  selCCEffFDHist =  new TH1D();
  selNCEffFDHist =  new TH1D();
  recoEffNDHist = new TH1D();
  recoCCEffFDHist = new TH1D();
  recoNCEffFDHist = new TH1D();
  extractNCFromCCFlux = new TH1D();
  fBeamMatrix =   new TH2D();
  
  //temporary histograms needed for true to reco 
  t2rCCnd = new TH2D("t2rCCnd","t2rCCnd", numBins, bins, numBins, bins );

  t2rCC = new TH2D("t2rCC","t2rCC",kNumTrueEnergyBins, &true_bins.at(0) , numBins, bins );
  t2rNC = new TH2D("t2rNC","t2rNC",kNumTrueEnergyBins, &true_bins.at(0) , numBins, bins );

  t2rCCRebin = new TH2D("t2rCCRebin","t2rCCRebin",numBins, bins , numBins, bins );
  t2rNCRebin = new TH2D("t2rNCRebin","t2rNCRebin",numBins, bins , numBins, bins );
  
  fSmoothWidth = 0; // Don't do any smoothing
  
  //needed for x section fitting
  bmFit = this ;
}
//----------------------------------------------------------------------



//......................................................................
NCExtrapolationBeamMatrix::~NCExtrapolationBeamMatrix(){}


const Registry& NCExtrapolationBeamMatrix::DefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  r.Set("FarNearSmoothingWidth", 0.00);
  r.Set("BeamMatrixFileLocation","/data/minos/pittam/beamdata/NewFluxHelpersRunIAllSkzpPiMinusNCUtilsBin.root");
  r.Set("SelectionEfficiencyFileLocation","/home/pittam/minos/devtest/EffMomCorrectionsWithNCNCUtilsBin.root");

  r.LockValues();
  return r;
}


//----------------------------------------------------------------------
//this method allows the user to pass in appropriate settings for the
//extrapolation such as fLocations of files, whether to generate the files
//again, etc.  the registry comes from the job module GetConfig method
void NCExtrapolationBeamMatrix::Prepare(const Registry& r)
{
  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "PREPARING BEAM MATRIX EXTRAPOLATION" << endl;

  // the base class takes care of turning the systematic parameters on/off
  // and setting their adjusted values if any

  //might need to change this to be passed in as part of the registry
  NCType::ECuts cutsType ;
  cutsType = NCType::kNCCuts ;
  SetFiducialMasses(cutsType) ;
    
  NCExtrapolation::Prepare(r);  

  int         tmpi;
  double      tmpd;
  const char* tmps;

  //  if(r.Get("UE3SqrVal",             tmpd)) fUE3Sqr                     = tmpd;

  if (r.Get("ExtractionType", tmpi)) fExtractionType=tmpi;

  // get location of beam matrix helpers
  if (r.Get("BeamMatrixFileLocation", tmps)) fBeamMatrixPath=tmps;
  GetBeamMatrix();
  // get selection and reconstruction efficiency histograms
  if (r.Get("SelectionEfficiencyFileLocation", tmps)) fSelEffFilePath=tmps;
  GetEfficiencyHistograms();
  if(r.Get("FarNearSmoothingWidth", tmpd)) fSmoothWidth = tmpd;

  //Logfile
  r.Get("FileName",tmps);
  string outFile(tmps);
  string::size_type start = outFile.rfind("/") + 1;
  outFile.erase(start,outFile.size());
  outFile += "ExtrapolationBeamMatrixFit.log"; 
  logFile = new ofstream(outFile.c_str(),ios::out | ios::app);
  
}


//----------------------------------------------------------------------
template<class T> void NCExtrapolationBeamMatrix::
HistSmoothHelper(T* h, bool is2D, double weight, double x, double y)
{
  if(fSmoothWidth == 0){
    // Shortcut all the below and just fill the histogram
    if(is2D)
      ((TH2D*)h)->Fill(x, y, weight);
    else
      h->Fill(x, weight);
    return;
  }

  // Paranoia over some scary text in the documentation for TAxis::GetBin
  h->SetBit(TH1::kCanRebin, false);

  TAxis* xaxis = h->GetXaxis();
  const int nbinsx = xaxis->GetNbins();

  // Warning - this number will only work in axis calls, not histogram calls
  const int xbin = xaxis->FindBin(x);
  const double lo = xaxis->GetBinLowEdge(xbin);
  const double hi = xaxis->GetBinUpEdge(xbin);

  const int ybin = is2D ? h->GetYaxis()->FindBin(y) : 0;

  const int histBin = xbin + (nbinsx+2)*ybin;

  const double dist_lo = x-lo;
  const double dist_hi = hi-x;

  // Don't smear events into the underflow bin
  const bool do_lo = xbin > 1 && dist_lo < fSmoothWidth*x;
  // Don't smear events into the overflow bin
  const bool do_hi = xbin < nbinsx && dist_hi < fSmoothWidth*x;

  double weight_here = 1;
  double weight_lo = 0;
  double weight_hi = 0;

  if(do_lo){
    const double transfer = .5*(1-dist_lo/(fSmoothWidth*x));
    weight_here -= transfer;
    weight_lo += transfer;
  }
  if(do_hi){
    const double transfer = .5*(1-dist_hi/(fSmoothWidth*x));
    weight_here -= transfer;
    weight_hi += transfer;
  }

  assert(weight_here >= .499);
  assert(weight_lo <= .501);
  assert(weight_hi <= .501);
  assert(weight_here <= 1);
  assert(weight_lo >= 0);
  assert(weight_hi >= 0);
  assert(fabs(weight_here+weight_lo+weight_hi-1) < .001);

  h->fArray[histBin] += weight*weight_here;

  if(do_lo){
    h->fArray[histBin-1] += weight*weight_lo;
  }
  if(do_hi){
    h->fArray[histBin+1] += weight*weight_hi;
  }
}

template void NCExtrapolationBeamMatrix::
HistSmoothHelper<TH1D>(TH1D*, bool, double, double, double);
template void NCExtrapolationBeamMatrix::
HistSmoothHelper<TH2D>(TH2D*, bool, double, double, double);

void NCExtrapolationBeamMatrix::FillHistogramSmoothed(TH1D* h, double x, double weight)
{
  HistSmoothHelper(h, false, weight, x);
}


void NCExtrapolationBeamMatrix::FillHistogramSmoothed2D(TH2D* h, double x, double y, double weight)
{
  HistSmoothHelper(h, true, weight, x, y);
}


//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::FillDataMCHistogramsNear(NCBeam* beam,
                                                         NCType::EEventType nccc,
                                                         int beamToUse)
{
  using Detector::kNear;

  //loop over NCEnergy bin objects in the near detector
  //to get the energies from the data events in each bin
  //and fill the fND_data histogram

  int selreco_ND = kNCSelRecoEnergy_ND;
  int reco_ND = kNCRecoEnergy_ND;
  int recoBkg = kNCRecoBkg_ND;
  int purity = kNCPurity_ND;
  int tau = kNCRecoEnergyTau_ND;
  int oscnue =  kNCRecoEnergyOscNue_ND ;
  int beamnue = kNCRecoEnergyBeamNueBkg_ND;

  int recototrue = kNCRecoTrueSig_ND;
  int recototrueBkg = kNCRecoTrueBkg_ND; 

  int recoData_ND = kNCDataRecoEnergy_ND;
  int trueE = kNCTrueEnergy_ND;
  int trueEBkg = kNCTrueEnergyBkg_ND;

  if(nccc == NCType::kCC){

    selreco_ND = kCCSelRecoEnergy_ND;
    reco_ND = kCCRecoEnergy_ND;
    recoBkg = kCCRecoBkg_ND;
    purity = kCCPurity_ND;
    tau = kCCRecoEnergyTau_ND;
    oscnue =  kCCRecoEnergyOscNue_ND ;
    beamnue = kCCRecoEnergyBeamNueBkg_ND;

    recototrue =  kCCRecoTrueSig_ND;
    recototrueBkg = kCCRecoTrueBkg_ND; 

    recoData_ND = kCCDataRecoEnergy_ND;
    trueE =  kCCTrueEnergy_ND;
    trueEBkg = kCCTrueEnergyBkg_ND;

  }
  
  //
  //Turns out that reseting is important
  //otherwise histograms pile up on each other
  //
  
  BMbybeams[beamToUse]->GetMCHists()[selreco_ND]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[reco_ND]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[recoBkg]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[purity]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[tau]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[beamnue]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[oscnue]->Reset( "ICE" );
  
  BMbybeams[beamToUse]->GetMCHists()[trueE]->Reset( "ICE" );
  BMbybeams[beamToUse]->GetMCHists()[trueEBkg]->Reset( "ICE" );
  
  BMbybeams[beamToUse]->Get2DHists()[recototrue]->Reset( "ICE" );
  BMbybeams[beamToUse]->Get2DHists()[recototrueBkg]->Reset( "ICE" );
  
  BMbybeams[beamToUse]->GetDataPredHists()[recoData_ND]->Reset( "ICE" );

  for(int b  = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
    NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

    const int dataSize       = bin->GetDataVectorSize();
    const int mcSize         = bin->GetMCSignalVectorSize();
    const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
    const int mcSize_tau     = bin->GetMCNuTauVectorSize();
    const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
    const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
        
    // The variables that get filled with event information
    double trueEnergy    = 0;
    double recoShowerE   = 0;
    double recoMuonE     = 0;
    double energy        = 0;
    double weight        = 0;
    int    flavor        = 0;

    for(int e = 0; e < dataSize; ++e){
      bin->GetDataInformation(energy, weight, e);
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetDataPredHists()[recoData_ND], energy, weight);
    }

    //------------------------------
    // numu
    //no nc weights * here as it doesn't matter only need CC spectrum to look right
    for(int e = 0; e < mcSize; ++e){
      bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed( BMbybeams[beamToUse]->GetMCHists()[reco_ND], recoShowerE + recoMuonE, weight);
      if(flavor == 14) FillHistogramSmoothed( BMbybeams[beamToUse]->GetMCHists()[purity], recoShowerE + recoMuonE, weight);    
      if(flavor == 14) FillHistogramSmoothed( BMbybeams[beamToUse]->GetMCHists()[trueE], trueEnergy, weight);    
      if(flavor == 14) FillHistogramSmoothed2D(BMbybeams[beamToUse]->Get2DHists()[recototrue], trueEnergy, recoShowerE + recoMuonE, weight);
    }

    for(int e = 0; e < mcSize_bkg; ++e){
      bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      if(nccc == NCType::kCC) weight *= BMbybeams[beamToUse]->GetmScaleNC().at(WhichFitParam(trueEnergy)) ;      
      // need to adjust CCBkg by the nd xsections... lets do this using only 1 2d histogram... meaning this gets made at the end of this method for the CC.
      // NC is fine as it is. Actually it doesn't even matter doesn't get used for anything.
      if(flavor == 14) FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[trueEBkg], trueEnergy, weight);    
      //need to add background to above histogram not subtract it. Including wrong sign in this one     
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[recoBkg], recoShowerE + recoMuonE, weight);             
      FillHistogramSmoothed2D(BMbybeams[beamToUse]->Get2DHists()[recototrueBkg], trueEnergy, recoShowerE + recoMuonE ,weight);
    }
  

    //------------------------------
    // nutau
    for(int e = 0; e < mcSize_tau; ++e){
      bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[tau], recoShowerE + recoMuonE, weight);
    }

    //------------------------------
    // beam nue
    for(int e = 0; e < mcSize_beamnue; ++e){
      bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[beamnue], recoShowerE + recoMuonE, weight);
    }

    //------------------------------
    // oscillated nue
    for(int e = 0; e < mcSize_oscnue; ++e){
      bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[oscnue], recoShowerE + recoMuonE, weight);
    }
  } // end for b
  
  
  // putting the data into the systematic beamTypes. Only needed for the LE beams. TODO change to correct runtype
  const NCBeam::Info beamInfoLEEEE = NCBeam::Info(BeamType::kL010z185i, NC::RunUtil::kRunAll);
  int nomBeamLE = fBMbybeams.find(beamInfoLEEEE)->second ;
  if(beam->GetBeamType() == BeamType::kL010z185i){
    for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
      BMbybeams[beamToUse]->GetDataPredHists()[recoData_ND]->SetBinContent(i,BMbybeams[nomBeamLE]->GetDataPredHists()[recoData_ND]->GetBinContent(i));
    }    
  }
  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    // purity correction histogram
    double denom = BMbybeams[beamToUse]->GetMCHists().at(reco_ND)->GetBinContent(i) + BMbybeams[beamToUse]->GetMCHists().at(recoBkg)->GetBinContent(i) + BMbybeams[beamToUse]->GetMCHists().at(beamnue)->GetBinContent(i);
    //make NC/cc selected histos
    BMbybeams[beamToUse]->GetMCHists().at(selreco_ND)->SetBinContent(i, denom ) ;
    if(denom > 0)
      BMbybeams[beamToUse]->GetMCHists().at(purity)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(purity)->GetBinContent(i) / denom ) ;
    else BMbybeams[beamToUse]->GetMCHists().at(purity)->SetBinContent(i,0);
  }
}

//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::FillDataMCHistogramsFar(NCBeam* beam,
                                                        NCType::EEventType nccc,
                                                        int beamToUse)
{
  using Detector::kFar;

  //loop over NCEnergy bin objects in the far detector
  double trueEnergy    = 0.;
  double recoShowerE   = 0.;
  double recoMuonE     = 0.;
  double weight        = 0.;
  int    flavor        = 0;

  int sig         = kNCSig;
  int bkg         = kNCBkg;
  int tau         = kNCTau;
  int beamnue     = kNCBeamNue;
  int oscnue      = kNCOscNue;
  int wrongsignbkg = kNCWrongSignBkg;

  int selreco_FD = kNCSelRecoEnergy_FD;
  int reco_FD = kNCRecoEnergy_FD;
  
  int trueE = kNCTrueEnergy_FD;
  int trueExsec = kNCTrueEnergyXsecFit_FD;
  int trueEBkg = kNCTrueEnergyBkg_FD;

  int recototrue = kNCTrueRecoSig_FD;
  int recototrueBkg = kNCTrueRecoBkg_FD; 
  
  int recoBkg = kNCRecoBkg_FD;
  int purity = kNCPurity_FD;

  int nue = kNCRecoEnergyBeamNueBkg_FD ;

  if(nccc == NCType::kCC){
    sig      = kCCSig;
    bkg      = kCCBkg;
    tau      = kCCTau;
    beamnue  = kCCBeamNue;
    oscnue   = kCCOscNue;
    wrongsignbkg = kCCWrongSignBkg;
  
    selreco_FD = kCCSelRecoEnergy_FD;
    reco_FD = kCCRecoEnergy_FD;
    
    trueE = kCCTrueEnergy_FD;
    trueExsec = kCCTrueEnergyXsecFit_FD;
    trueEBkg = kCCTrueEnergyBkg_FD;

    recototrue = kCCTrueRecoSig_FD;
    recototrueBkg = kCCTrueRecoBkg_FD; 

    recoBkg = kCCRecoBkg_FD;
    purity = kCCPurity_FD;

    nue = kCCRecoEnergyBeamNueBkg_FD ;
}

   BMbybeams[beamToUse]->Get2DOscHists()[sig]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DOscHists()[bkg]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DOscHists()[tau]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DOscHists()[beamnue]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DOscHists()[oscnue]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DOscHists()[wrongsignbkg]->Reset( "ICE" );

   BMbybeams[beamToUse]->GetMCHists()[selreco_FD]->Reset( "ICE" );
   BMbybeams[beamToUse]->GetMCHists()[reco_FD]->Reset( "ICE" );
   
   BMbybeams[beamToUse]->GetMCHists()[trueE]->Reset( "ICE" );
   BMbybeams[beamToUse]->GetMCHists()[trueExsec]->Reset( "ICE" );
   BMbybeams[beamToUse]->GetMCHists()[trueEBkg]->Reset( "ICE" );
   
   BMbybeams[beamToUse]->Get2DHists()[recototrue]->Reset( "ICE" );
   BMbybeams[beamToUse]->Get2DHists()[recototrueBkg]->Reset( "ICE" );

   BMbybeams[beamToUse]->GetMCHists()[recoBkg]->Reset( "ICE" );
   BMbybeams[beamToUse]->GetMCHists()[purity]->Reset( "ICE" );
   
   BMbybeams[beamToUse]->GetMCHists()[nue]->Reset( "ICE" );

   const int B = beam->GetNumberEnergyBins(nccc);
   for(int b = 0; b < B; ++b){
     NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

     const int mcSize         = bin->GetMCSignalVectorSize();
     const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
     const int mcSize_tau     = bin->GetMCNuTauVectorSize();
     const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
     const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
     
     //numu
     for(int e = 0; e < mcSize; ++e){
       bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
       if(nccc == NCType::kNC) recoMuonE = 0;
       if(flavor == 14) FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[trueE], trueEnergy, weight);           
       //The weights for the xsection go into the true to reco histogram and are applied at the true to reco stage. 
       //don't need to do the other 2d ones as they are done on the nominal mc. The purity hist is only used in the unosc prediction at the moment. Will need to be adjusted if not though
       //Relevant Backgrounds are adjusted for xsecs in the DoOscillations method
       //When fitting for a NC Clean systematic it only shifts the nc selected and wrong sign bkg NOT cc bkg... small effect.
             
      if(flavor == 14)FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[sig], recoShowerE + recoMuonE, trueEnergy, weight);      
      if(flavor != 14)FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[wrongsignbkg], recoShowerE + recoMuonE, trueEnergy, weight);
      if(nccc == NCType::kNC) weight *= BMbybeams[beamToUse]->GetmScaleNC().at(WhichFitParam(trueEnergy));
      if(flavor == 14) FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[purity], recoShowerE + recoMuonE, weight);    
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[reco_FD], recoShowerE + recoMuonE, weight);
      if(flavor == 14) FillHistogramSmoothed2D(BMbybeams[beamToUse]->Get2DHists()[recototrue], trueEnergy, recoShowerE + recoMuonE, weight);
      if(flavor == 14) FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[trueExsec], trueEnergy, weight);           
      

    }

    for(int e = 0; e < mcSize_bkg; ++e){
      bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      if(flavor == 14) FillHistogramSmoothed2D(BMbybeams[beamToUse]->Get2DHists()[recototrueBkg], trueEnergy, recoShowerE + recoMuonE, weight);
      if(flavor == 14) FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[trueEBkg], trueEnergy, weight);     
      if(nccc == NCType::kCC) weight *= BMbybeams[beamToUse]->GetmScaleNC().at(WhichFitParam(trueEnergy));
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[recoBkg], recoShowerE + recoMuonE, weight);      
      //Weighting this here rather than DoOsc method
      FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[bkg], recoShowerE + recoMuonE, trueEnergy, weight);
    }
    
    //nutau
    for(int e = 0; e < mcSize_tau; ++e){
      bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[tau], recoShowerE + recoMuonE, trueEnergy, weight);
    }

    //beamnue
    for(int e = 0; e < mcSize_beamnue; ++e){
      bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetMCHists()[nue], recoShowerE + recoMuonE, weight);
      FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[beamnue], recoShowerE + recoMuonE, trueEnergy, weight);

    }
    
    //oscnue
    for(int e = 0; e < mcSize_oscnue; ++e){
      bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D( BMbybeams[beamToUse]->Get2DOscHists()[oscnue], recoShowerE + recoMuonE, trueEnergy, weight);
    }
    
    const int N = bin->GetDataVectorSize();
    for(int n = 0; n < N; ++n){
      double energy, weight;
      bin->GetDataInformation(energy, weight, n);
      FillHistogramSmoothed(BMbybeams[beamToUse]->GetFDDataHists()[sig], energy, weight);
    } // end for n
  } // end for b

  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    double denom = BMbybeams[beamToUse]->GetMCHists().at(reco_FD)->GetBinContent(i) + BMbybeams[beamToUse]->GetMCHists().at(recoBkg)->GetBinContent(i) + BMbybeams[beamToUse]->GetMCHists().at(nue)->GetBinContent(i) ;
    BMbybeams[beamToUse]->GetMCHists().at(selreco_FD)->SetBinContent(i, denom ) ;    
    if(denom > 0)
      BMbybeams[beamToUse]->GetMCHists().at(purity)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(purity)->GetBinContent(i) / denom ) ;
    else BMbybeams[beamToUse]->GetMCHists().at(purity)->SetBinContent(i,0);
  }


}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::
ConstructFarSpectrum(const NC::OscProb::OscPars* pars, int beamsToUse)
{
  DoPurityCorrectionND(beamsToUse);
  DoRecoToTrueND(beamsToUse);
  DoEfficencyCorrectionND(beamsToUse);
  DoXsectionCorrectionND(beamsToUse);
  DoBeamMatrixExtrapolation(beamsToUse);
  GetNCFromCCFlux(beamsToUse);
  DoXsectionCorrectionFD(beamsToUse);
  DoEfficencyCorrectionFD(beamsToUse);
  DoUnoscTrueToRecoFD(beamsToUse);
  DoUnoscPurityCorrectionFD(beamsToUse);
  //need to call unosc before oscillations
  DoOscillations(pars, beamsToUse);
}
//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::
FindSpectraForPars(const NC::OscProb::OscPars* oscPars,
                   const NC::SystPars& systPars_org,
                   vector<NCBeam::Info> beamsToUse,
                   vector<TH1*>& expvec,
                   vector<TH1*>& obsvec)
{


  //not sure where to put this...
  //has to be before first call to filldatamchistograms
  //but after addevent loops

  // these are taken care of by clever weighting of NCBeams
  fFDPOTData = 1.0 ;
  fFDPOTMC = 1.0 ;
  
  fNDPOTData = 1.0 ;
  fNDPOTMC =  1.0 ;

  int icount = 0;  
  int ibm = 0;

  static bool doneFillHists = false ;
  if(!doneFillHists) {

    for (Int_t i=0; i<4; i++) {
      m0ScaleNC.push_back(1.0);
      m0ScaleNCErr.push_back(0.005);

      mScaleNC.push_back(1.0);
      mScaleNCErr.push_back(0.005);

    }
    //this is for use in the nd xsection fitting
    for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
      const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
      const NCBeam::Info beamInfo = key.first;
      const Detector::Detector_t det = key.second;    
      NCBeam* beam = it->second;        
      if(det == Detector::kNear) beamsNear.push_back(beam); 
      MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Beam Description " << beamInfo.GetDescription() << " det " <<  beam->GetDetector() <<endl;

    }

    //filling the vector of BeamMatrix Objects ibm is index for Map
    for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
      
      const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
      const NCBeam::Info beamInfo = key.first;
      const Detector::Detector_t det = key.second;

      NCBeam* beam = it->second;        
      
      string name = (string)beam->GetInfo().GetDescription();
      string title = "title_"+name;
      string bmname = "beamMatrix_"+name;
      
      //only 1 BeamMatrix object for ND and FD
      if(icount%2 == 0){ 
        BMbybeams.push_back(new BeamMatrix(bmname,title,name));
        fBMbybeams[beamInfo] = ibm;
      }
      
      if( beam->GetDetector() == Detector::kNear){
        //xsection fit is done here before FillDataMCHistograms
        if(beam->GetBeamType() == BeamType::kL010z185i) doNDXsectionFit(1, beam, ibm);
        if(fUseNC) FillDataMCHistogramsNear(beam, NCType::kNC, ibm);
        if(fUseCC) FillDataMCHistogramsNear(beam, NCType::kCC, ibm);
        MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Filled Beam with Description " << name << " det  " << det << " ibm " << ibm 
                                                     << " BeamsNear.at(ibm)->GetDescription() " << beamsNear.at(ibm)->GetInfo().GetDescription() << endl;
      }
        
      if( beam->GetDetector() == Detector::kFar){
        BMbybeams[ibm]->GetFDDataHists().at(kNCSig)->Reset("ICE");
        BMbybeams[ibm]->GetFDDataHists().at(kCCSig)->Reset("ICE");
        if(fUseNC) FillDataMCHistogramsFar(beam, NCType::kNC, ibm);
        if(fUseCC) FillDataMCHistogramsFar(beam, NCType::kCC, ibm);
        
        MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Filled Beam with Description " << name << " det  " << det << " ibm " << ibm <<endl ;     
      }
      
      if(icount%2 == 1) ibm++; 
      icount++ ;
    }
  
  }
  //only need to do nd xsec fit and fill hists once
  doneFillHists = true ;

  // making fd prediciton for beams asked for
  for (int i=0; i<(int)beamsToUse.size(); ++i) {
    
    NCBeam::Info beamInfo = beamsToUse[i];
    
    //getting beams and index for BeamMatrix object that correpsonds to those beams
    int bmeq = fBMbybeams.find(beamInfo)->second;
    
    NCBeam* beam_fd = GetBeam(Detector::kFar , beamInfo);
    NCBeam* beam_nd = GetBeam(Detector::kNear, beamInfo);

    assert(beam_fd);
    assert(beam_nd);
    
    //only want to extrapolate LE beams. Others are only used for ND fitting
    if(beam_fd->GetBeamType() !=  BeamType::kL010z185i) continue ;

    if(fDoSystematicInterpolation){
      InitializeInterpolator(oscPars);
    }
    
    if(!fInterpolator){
      
      BMbybeams[bmeq]->GetFitHists().at(kNCSig)->Reset("ICE"); //integral contents error
      BMbybeams[bmeq]->GetFitHists().at(kCCSig)->Reset("ICE");
      
      ConstructFarSpectrum(oscPars,bmeq);
      
      if(fUseNC){
        expvec.push_back(BMbybeams[bmeq]->GetFitHists()[kNCSig]);
        obsvec.push_back(BMbybeams[bmeq]->GetFDDataHists()[kNCSig]);
      }
      
      if(fUseCC){
        expvec.push_back(BMbybeams[bmeq]->GetFitHists()[kCCSig]);
        obsvec.push_back(BMbybeams[bmeq]->GetFDDataHists()[kCCSig]);
      }      
      return ;
    }
   
    const vector<double> shift = fCoordConv.VectorFromSystPars(systPars_org);
    vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);
    
    BMbybeams[i]->GetFitHists().at(kNCSig)->Reset("ICE"); //integral contents error
    BMbybeams[i]->GetFitHists().at(kCCSig)->Reset("ICE");

    ConstructFarSpectrum(oscPars ,bmeq);    
    
    MultiplyFast(BMbybeams[bmeq]->GetFitHists()[kNCSig], interp[0]);
    MultiplyFast(BMbybeams[bmeq]->GetFitHists()[kCCSig], interp[2]);
    
    if(fUseNC){
      expvec.push_back(BMbybeams[bmeq]->GetFitHists()[kNCSig]);
      obsvec.push_back(BMbybeams[bmeq]->GetFDDataHists()[kNCSig]);
    }
    
    if(fUseCC){
      expvec.push_back(BMbybeams[bmeq]->GetFitHists()[kCCSig]);
      obsvec.push_back(BMbybeams[bmeq]->GetFDDataHists()[kCCSig]);
    }
  }
  ++fNCalls;
  return;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::CleanupSpectra(vector<TH1*> /*exps*/,
                                            vector<TH1*> /*obss*/)
{
  // No need to clear anything. We don't allocate anything new each time
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::
WriteResources(const NC::OscProb::OscPars* trueOscPars)
{
  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "nCalls: " << fNCalls << endl;

  // only writing out RunI at the moment
  NCBeam::Info bmInfo  = NCBeam::Info(BeamType::kL010z185i,NC::RunUtil::kRunI) ;
  int bmequiv = fBMbybeams.find(bmInfo)->second;
  
  NC::SystPars systPars = GetBestFitSysts();
  const NC::OscProb::OscPars* oscPars = GetBestFitOscPars();
  
  if(fDoSystematicInterpolation){
    const vector<double> shift = fCoordConv.VectorFromSystPars(systPars);
    vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);
    
    ConstructFarSpectrum(oscPars, bmequiv);    
    
    MultiplyFast(BMbybeams[bmequiv]->GetFitHists()[kNCSig], interp[0]);
    MultiplyFast(BMbybeams[bmequiv]->GetFitHists()[kCCSig], interp[2]);
    
  }
  
  else ConstructFarSpectrum(oscPars, bmequiv);

  delete oscPars;
    
  if(fUseNC) FillResultHistograms(bmInfo, NCType::kNC, fNominalBM, fOscBM);
  if(fUseCC) FillResultHistograms(bmInfo, NCType::kCC, fNominalBM, fOscBM);  

  NCExtrapolation::WriteResources(trueOscPars);
  // get a pointer to the current directory
  // this is one of the output files

  TDirectory* saveDir = gDirectory;

  //writing out histograms. Any others wanted can be easily added ME or HE beam for example
  for(unsigned int i = 0; i < 2; ++i) BMbybeams[bmequiv]->GetFDDataHists()[i]->Write();
  for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetMCHists().size(); ++i) BMbybeams[bmequiv]->GetMCHists().at(i)->Write();
  for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetDataPredHists().size(); ++i) BMbybeams[bmequiv]->GetDataPredHists().at(i)->Write();
  for(unsigned int i = 0; i < BMbybeams[bmequiv]->Get2DHists().size(); ++i)  BMbybeams[bmequiv]->Get2DHists().at(i)->Write();
  for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetFitHists().size();  ++i)  BMbybeams[bmequiv]->GetFitHists().at(i)->Write();
  for(unsigned int i = 0; i < ND_XSectionFitRew.size();  ++i)  ND_XSectionFitRew[i]->Write();
  for(unsigned int i = 0; i < ND_XSectionFit.size();  ++i)  ND_XSectionFit[i]->Write();
  
  selEffNDHist->Write() ;
  selCCEffFDHist->Write() ;
  selNCEffFDHist->Write() ;
  recoEffNDHist->Write() ;
  recoCCEffFDHist->Write() ;
  recoNCEffFDHist->Write() ;
  extractNCFromCCFlux->Write() ;
  saveDir->cd();

  //writing the result of the x-section fit written into the file
  TString parname ;
  TString parout ;
  gDirectory->mkdir("xSecFit", "xSecFit")->cd();
  Registry* xsecFit = new Registry;
  xsecFit->SetName("xSecFit_registry");
  xsecFit->UnLockKeys();
  xsecFit->UnLockValues();
  for (Int_t i=0; i<4; i++){ 
    parname = (TString) Form("Par %i: ",i);
    parout = (TString) Form("%e +%e %e +/- %e, Corr: %e",mScaleNC[i],mScalePos[i],mScaleNeg[i],mScaleNCErr[i],corr[i]);
    xsecFit->Set(parname,parout);
  }
  xsecFit->Write();
  saveDir->cd();
  delete xsecFit;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------

void NCExtrapolationBeamMatrix::
FillResultHistograms(NCBeam::Info beamInfo,
                     NCType::EEventType nccc,
                     vector<vector<double> >& nominal,
                     vector<vector<double> >& osc){
  
  NCBeam* beam = GetBeam(Detector::kFar, beamInfo);
  int bmequiv = fBMbybeams.find(beamInfo)->second;

  int sig     = kNCSig;
  int bkg     = kNCBkg;
  int tau     = kNCTau;
  int beamnue = kNCBeamNue;
  int oscnue  = kNCOscNue;
  int wrongsignbkg = kNCWrongSignBkg;

  if(nccc == NCType::kCC){
    sig     = kCCSig;
    bkg     = kCCBkg;
    tau     = kCCTau;
    beamnue = kCCBeamNue;
    oscnue  = kCCOscNue;
    wrongsignbkg = kCCWrongSignBkg;
  }

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

    double eBeamVal = nominal[beamnue][i];
    double nooscVal = nominal[sig][i] + nominal[bkg][i] + nominal[wrongsignbkg][i] + eBeamVal;
    double fitVal   = BMbybeams[bmequiv]->GetFitHists()[sig]->GetBinContent(i+1);
    double bkgVal   = osc[bkg][i] ;
    double tauVal   = osc[tau][i];
    double eOscVal  = osc[oscnue][i]  ;

    if(nooscVal > 0 )
      beam->GetMCHistogram(nccc)->SetBinContent(i+1, nooscVal);
    else beam->GetMCHistogram(nccc)->SetBinContent(i+1, 0);
    if(fitVal > 0)
      beam->GetMCFitHistogram(nccc)->SetBinContent(i+1, fitVal);
    else beam->GetMCFitHistogram(nccc)->SetBinContent(i+1, 0);
    if(bkgVal > 0){
      beam->GetMCBackgroundHistogram(nccc)->SetBinContent(i+1, bkgVal);
    }
    else  beam->GetMCBackgroundHistogram(nccc)->SetBinContent(i+1, 0);
    if(tauVal > 0)
      beam->GetMCFitNuMuToNuTauHistogram(nccc)->SetBinContent(i+1, tauVal);
    else beam->GetMCFitNuMuToNuTauHistogram(nccc)->SetBinContent(i+1,0);
    if(eBeamVal > 0)
      beam->GetMCFitNuEToNuEHistogram(nccc)->SetBinContent(i+1, eBeamVal);
    else beam->GetMCFitNuEToNuEHistogram(nccc)->SetBinContent(i+1, 0);
    if( eOscVal > 0)
      beam->GetMCFitNuMuToNuEHistogram(nccc)->SetBinContent(i+1, eOscVal);
    else  beam->GetMCFitNuMuToNuEHistogram(nccc)->SetBinContent(i+1, 0);

    double bkgValnofit   = nominal[bkg][i];
    double tauValnofit   = nominal[tau][i];
    double eOscValnofit  = nominal[oscnue][i];
    double eBeamValnofit = osc[beamnue][i];

    if(bkgValnofit > 0)
      beam->GetMCNoFitBackgroundHistogram(nccc)->SetBinContent(i+1,bkgValnofit);
    else beam->GetMCNoFitBackgroundHistogram(nccc)->SetBinContent(i+1,0);
    
    if( tauValnofit > 0)
      beam->GetMCNoFitNuMuToNuTauHistogram(nccc)->SetBinContent(i+1,tauValnofit);
    else  beam->GetMCNoFitNuMuToNuTauHistogram(nccc)->SetBinContent(i+1,0);
    
    if(eOscValnofit > 0)
      beam->GetMCNoFitNuMuToNuEHistogram(nccc)->SetBinContent(i+1,tauValnofit);
    else  beam->GetMCNoFitNuMuToNuEHistogram(nccc)->SetBinContent(i+1,0);
    
    if(eBeamValnofit > 0)
      beam->GetMCNoFitNuEToNuEHistogram(nccc)->SetBinContent(i+1,eBeamValnofit);
    else  beam->GetMCNoFitNuEToNuEHistogram(nccc)->SetBinContent(i+1,0);
    
  }

  // We filled the histograms manually, and don't need NCBeam's automatic stuff
  // to happen. So inform it of that fact.
  beam->MarkHistogramsFilled();
}
//----------------------------------------------------------------------

//----------------------------------------------------------------------
bool NCExtrapolationBeamMatrix::EnableNearToFarExtrapolation(bool enable)
{
  const bool ret = fDoExtrapolation;
  fDoExtrapolation = enable;
  return ret;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// This interface is most convenient for internal use
// I don't know what I need this for...
// It certainly doesn't return the NearSpectra
vector<const TH1*> NCExtrapolationBeamMatrix::
GetNearMCSpectra(vector<NCBeam::Info> beamsToUse)
{
  
  vector<const TH1*> ret;
  for (int i=0; i<(int)beamsToUse.size(); ++i) {
    NCBeam::Info beamInfo = beamsToUse[i];
    
    int bmequiv = fBMbybeams.find(beamInfo)->second;
    NCBeam* beam_nd= GetBeam(Detector::kNear, beamInfo);
    
    if(beam_nd->GetBeamType() !=  BeamType::kL010z185i) continue ;

    // I suspect that what you return here makes no difference. 
    //It does have to be return the correct number of histograms though otherwise seg faults
    
    if(fUseNC) ret.push_back(BMbybeams[bmequiv]->GetMCHists().at(kNCSelRecoEnergy_ND));
    if(fUseCC) ret.push_back(BMbybeams[bmequiv]->GetMCHists().at(kCCSelRecoEnergy_ND));
   
  }
  return ret;
  
}
//----------------------------------------------------------------------

//----------------------------------------------------------------------
void NCExtrapolationBeamMatrix::
DoPurityCorrectionND(int beamToUse){
 
  for(int i = 0; i < kNumEnergyBinsFar+1; ++i ){
    BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCSelRecoEnergy_ND)->GetBinContent(i)
                                                                                  * BMbybeams[beamToUse]->GetMCHists().at(kCCPurity_ND)->GetBinContent(i));
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_ND)->SetBinContent(i, BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergy_ND)->GetBinContent(i) 
                                                                                            * BMbybeams[beamToUse]->GetMCHists().at(kCCPurity_ND)->GetBinContent(i));

  }
  return;  
}

void NCExtrapolationBeamMatrix::
DoRecoToTrueND(int beamToUse){

  vector<double> tempMC(kNumEnergyBinsFar+1,0) ;
  vector<double> tempData(kNumEnergyBinsFar+1,0) ;
  
  t2rCCnd = (TH2D*)CloneFast(BMbybeams[beamToUse]->Get2DHists()[kCCRecoTrueSig_ND]);

  double bincontent = 0;  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//true
      // I think I want to use the ndmc to normalise to get cancellations 
      // Normalise reco to true matrix
      if( BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i) > 0){
        bincontent = (t2rCCnd->GetBinContent(j,i) / BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i));
        t2rCCnd->SetBinContent(j,i, bincontent);
      }
      else t2rCCnd->SetBinContent(j,i,0);
    }
  }
  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//true
      // correct reco to true
      tempMC[j-1] += (BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i) * t2rCCnd->GetBinContent(j,i) );
      tempData[j-1] += (BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_ND)->GetBinContent(i) * t2rCCnd->GetBinContent(j,i) );
    }
  }
  
  for(int i=1;i<=kNumEnergyBinsFar;i++){
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyPurCorr_ND)->SetBinContent(i,tempMC[i-1]);
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyPurCorr_ND)->SetBinContent(i,tempData[i-1]);
    
  }

  delete t2rCCnd ;
  tempMC.clear() ;
  tempData.clear() ;
  return;  
}

void NCExtrapolationBeamMatrix::
GetEfficiencyHistograms(){

  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Getting Efficiency Histograms from: " << fSelEffFilePath << " Extraction Type " << fExtractionType << endl;
  TFile* selEffFile = new TFile(fSelEffFilePath,"READ");
  if(selEffFile){
    
    extractNCFromCCFlux = (TH1D*)selEffFile->Get("extractNCFromCCFlux");
    
    recoEffNDHist = (TH1D*)selEffFile->Get("recoEffHistND");
    recoCCEffFDHist = (TH1D*)selEffFile->Get("recoCCEffHistFD");
    recoNCEffFDHist = (TH1D*)selEffFile->Get("recoNCEffHistFD");

    switch( fExtractionType){
    case NCType::kNCCCExtractionCuts:
      selEffNDHist = (TH1D*)selEffFile->Get("selCorrectionTOmND");
      selCCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionCCTOmFD"); 
      selNCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionNCTOmFD"); 
      break;
    case NCType::kNCCCExtractionANN:
    case NCType::kNCCCExtractionANNNear:
    case NCType::kNCCCExtractionANNFar:
      selEffNDHist = (TH1D*)selEffFile->Get("selCorrectionRPannND"); 
      selCCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionCCRPannFD"); 
      selNCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionNCRPannFD"); 
      break;
    case NCType::kNCCCExtractionkNN:
      selEffNDHist = (TH1D*)selEffFile->Get("selCorrectionROND"); 
      selCCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionCCROFD"); 
      selNCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionNCROFD"); 
      break;
    case NCType::kNCCCExtractionMDA:
      selEffNDHist = (TH1D*)selEffFile->Get("selCorrectionASND"); 
      selCCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionCCASFD"); 
      selNCEffFDHist = (TH1D*)selEffFile->Get("selCorrectionNCASFD"); 
      break;
    default:
      selEffNDHist = 0 ;
      selCCEffFDHist = 0;
      selNCEffFDHist = 0;
    }    
  }
  
  else{ 
    MSG("NCExtrapolationBeamMatrix", Msg::kError) <<  "Can't get Selection Efficiency Histograms from " << fSelEffFilePath  <<  " Extraction Type " <<  fExtractionType << endl;
  }
  
}

void NCExtrapolationBeamMatrix::
DoEfficencyCorrectionND(int beamToUse){
  
  //do selection then reconstruction in the ND
     
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    
    // selection efficiency correction

    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_ND)->SetBinContent(i,selEffNDHist->GetBinContent(i));
        
    if(BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_ND)->GetBinContent(i) > 0){
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyPurCorr_ND)->GetBinContent(i) / BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_ND)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyPurCorr_ND)->GetBinContent(i) / BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_ND)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->SetBinContent(i,0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->SetBinContent(i,0);
    }


  // reconstruction efficiency correction    
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_ND)->SetBinContent(i,recoEffNDHist->GetBinContent(i));
    
    if(BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_ND)->GetBinContent(i) > 0){
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->GetBinContent(i) / BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_ND)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->GetBinContent(i) / BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_ND)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->SetBinContent(i,0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->SetBinContent(i,0);
    }
    

  }
}
void NCExtrapolationBeamMatrix::
DoXsectionCorrectionND(int beamToUse){
  
  // Doesn't actually do Xsections as they are folded into the BeamMatrix
  // it does do the fiducial mass and pot correction

  MAXMSG("NCExtrapolationBeamMatrix", Msg::kInfo, 1) << "ND Fid Vol " << fNDFidVolMass 
                                               << " ND POT MC " << fNDPOTMC
                                               << " ND POT Data " << fNDPOTData
                                               << endl;

  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyFlux_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i)*fFluxPOT/(fNDFidVolMass*fNDPOTMC));
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyFlux_ND)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->GetBinContent(i)/(fNDFidVolMass*fNDPOTData));
  }
}

void NCExtrapolationBeamMatrix::
GetBeamMatrix(){

  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "Getting Beam Matrix from: " << fBeamMatrixPath << endl;
  TFile* beamMatrixFile = new TFile(fBeamMatrixPath,"READ");
  if(beamMatrixFile){
    fBeamMatrix =  (TH2D*)beamMatrixFile->Get("FDVsNDMatrixXSecRW");
  }
  
  else{ 
    MSG("NCExtrapolationBeamMatrix", Msg::kError) <<  "Can't get BeamMatrix from " << fBeamMatrixPath  << endl;
  }
}

void NCExtrapolationBeamMatrix::
DoBeamMatrixExtrapolation(int beamToUse){

  int beamMatrixNearBins = fBeamMatrix->GetXaxis()->GetNbins();
  int beamMatrixFarBins = fBeamMatrix->GetYaxis()->GetNbins();

  vector<double> temp(beamMatrixNearBins+1,0) ;
  vector<double> tempData(beamMatrixNearBins+1,0) ;
  
  for(int i = 1; i < beamMatrixNearBins+1; ++i ){
    for(int j = 1; j <beamMatrixFarBins+1 ; ++j ){      
      temp[j] += (BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyFlux_ND)->GetBinContent(i) * fBeamMatrix->GetCellContent( i,j ));
      tempData[j] += (BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyFlux_ND)->GetBinContent(i) * fBeamMatrix->GetCellContent( i,j ));
    }
  }

  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyMatrix_FD)->SetBinContent(i,temp[i]);
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->SetBinContent(i,tempData[i]);

    
  }  
  temp.clear() ;
  tempData.clear() ;
}

void NCExtrapolationBeamMatrix::
GetNCFromCCFlux(int beamToUse){
  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){

    //not really sure where to put this either here or any of the following steps.
    //But the histograms already exist so it may as well be here. 
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyMatrix_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyMatrix_FD)->GetBinContent(i) 
                                                                                 * extractNCFromCCFlux->GetBinContent(i));
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyMatrix_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->GetBinContent(i) 
                                                                                           * extractNCFromCCFlux->GetBinContent(i));    
  }
}
void NCExtrapolationBeamMatrix::
SetFiducialMasses(NCType::ECuts cutsType){
  
  fFluxPOT = 1.0;
  
  switch(cutsType){
  case NCType::kCCCuts:
    fNDFidVolMass = -9999.9 ;
    fFDFidVolMass = -9999.9 ;
    break;
  case NCType::kNCCuts:
    fNDFidVolMass = 26874.0; //kg
    fFDFidVolMass = 3726961.; //kg
    break;
  case NCType::kNCCCFidCuts:
    fNDFidVolMass = -9999.9 ;
    fFDFidVolMass = -9999.9 ;   
    break;
  default:
    fNDFidVolMass = -9999.9 ;
    fFDFidVolMass = -9999.9 ;   
  }
}

void NCExtrapolationBeamMatrix::
DoXsectionCorrectionFD(int beamToUse){

  // Doesn't actually do Xsections as they are folded into the BeamMatrix
  // it does do the fiducial mass and pot correction
  
  MAXMSG("NCExtrapolationBeamMatrix", Msg::kInfo, 1) << "FD Fid Vol " << fFDFidVolMass 
                                                      << " FD POT MC " << fFDPOTMC
                                                      << " FD POT Data " << fFDPOTData
                                                      << endl;  
  

  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    
    //CC
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyFlux_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyMatrix_FD)->GetBinContent(i)
                                                                               *(fFDFidVolMass*fFDPOTMC)/fFluxPOT);
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyFlux_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->GetBinContent(i)
                                                                                         *(fFDFidVolMass*fFDPOTData));
    
    //NC
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyFlux_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyMatrix_FD)->GetBinContent(i)
                                                                               *(fFDFidVolMass*fFDPOTMC)/fFluxPOT);
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyFlux_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyMatrix_FD)->GetBinContent(i)
                                                                                         *(fFDFidVolMass*fFDPOTData));
    
  }
}

void NCExtrapolationBeamMatrix::
DoEfficencyCorrectionFD(int beamToUse){
 
  //
  // Do reconstruction then selection efficiency in the Far Detector
  //
  
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
    
   
    // reconstruction efficiency correction    
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_FD)->SetBinContent(i,recoCCEffFDHist->GetBinContent(i));
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueRecoEfficiency_FD)->SetBinContent(i,recoNCEffFDHist->GetBinContent(i));    

    // CC correction
    if(BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_FD)->GetBinContent(i) > 0){
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyFlux_FD)->GetBinContent(i) 
                                                                                    * BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_FD)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyFlux_FD)->GetBinContent(i) 
                                                                                              * BMbybeams[beamToUse]->GetMCHists().at(kCCTrueRecoEfficiency_FD)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->SetBinContent(i,0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->SetBinContent(i,0);
    }
    
    // NC correction
    if(BMbybeams[beamToUse]->GetMCHists().at(kNCTrueRecoEfficiency_FD)->GetBinContent(i) > 0){
      BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyFlux_FD)->GetBinContent(i) 
                                                                                    * BMbybeams[beamToUse]->GetMCHists().at(kNCTrueRecoEfficiency_FD)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyFlux_FD)->GetBinContent(i) 
                                                                                              * BMbybeams[beamToUse]->GetMCHists().at(kNCTrueRecoEfficiency_FD)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->SetBinContent(i,0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->SetBinContent(i,0);
    }

    
    // selection efficiency correction
    BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_FD)->SetBinContent(i,selCCEffFDHist->GetBinContent(i));
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueSelEfficiency_FD)->SetBinContent(i,selNCEffFDHist->GetBinContent(i));
    
    //CC
    if(BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_FD)->GetBinContent(i) > 0){
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->GetBinContent(i) 
                                                                                       * BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_FD)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->GetBinContent(i) 
                                                                                                 * BMbybeams[beamToUse]->GetMCHists().at(kCCTrueSelEfficiency_FD)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
    }
  
    //NC
    if(BMbybeams[beamToUse]->GetMCHists().at(kNCTrueSelEfficiency_FD)->GetBinContent(i) > 0){
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->GetBinContent(i) 
                                                                                     * BMbybeams[beamToUse]->GetMCHists().at(kNCTrueSelEfficiency_FD)->GetBinContent(i) );
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->GetBinContent(i) 
                                                                                               * BMbybeams[beamToUse]->GetMCHists().at(kNCTrueSelEfficiency_FD)->GetBinContent(i) );
    }
    else{
    BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
    }

  }  
}

void NCExtrapolationBeamMatrix::
DoUnoscTrueToRecoFD(int beamToUse){
  
  
  vector<double> tempCCMC(kNumEnergyBinsFar+1,0) ;
  vector<double> tempNCMC(kNumEnergyBinsFar+1,0) ;
  vector<double> tempCCData(kNumEnergyBinsFar+1,0) ;
  vector<double> tempNCData(kNumEnergyBinsFar+1,0) ;

  //need to use nominal for normalisation or Systematics cancel at true to reco stage
  const NCBeam::Info beamInfoLE = NCBeam::Info(BeamType::kL010z185i, NC::RunUtil::kRunI);
  int leNomBeam = fBMbybeams.find(beamInfoLE)->second ;

  //should be ok to use these here
  t2rCCRebin->Reset("ICE");
  t2rNCRebin->Reset("ICE");
  
  t2rCC = (TH2D*)CloneFast(BMbybeams[beamToUse]->Get2DHists()[kCCTrueRecoSig_FD]);
  t2rNC = (TH2D*)CloneFast(BMbybeams[beamToUse]->Get2DHists()[kNCTrueRecoSig_FD]);
  
  //see DoOscillations for the 4 step plan
  //step 1
  //unosc so no need for this step
  //step 2

  rebinRecoTrueToRecoRecoHist(t2rCC, t2rCCRebin); 
  rebinRecoTrueToRecoRecoHist(t2rNC, t2rNCRebin); 

  //step 3
  normaliseRecoToTrue(t2rCCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kCCTrueEnergy_FD ) );
  normaliseRecoToTrue(t2rNCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kNCTrueEnergy_FD ) );

  //step 4
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//true
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
      // correct reco to true CC
      tempCCMC[j-1] += (BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j) );
      tempCCData[j-1] += (BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j));
      
      // correct reco to true NC
      tempNCMC[j-1] += (BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j) );
      tempNCData[j-1] += (BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j) );
    }
  }
  
  for(int i = 0; i < kNumEnergyBinsFar; ++i ){//reco
    BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempCCMC[i]);
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempCCData[i]);
    
    // correct reco to true NC
    BMbybeams[beamToUse]->GetMCHists().at(kNCRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempNCMC[i]);
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempNCData[i]);
  }

  tempCCMC.clear() ;
  tempCCData.clear() ;
  tempNCMC.clear() ;
  tempNCData.clear() ;
  
  delete t2rCC;
  delete t2rNC;
  
  return;  
}

void NCExtrapolationBeamMatrix::
DoUnoscPurityCorrectionFD(int beamToUse){
 
  for(int i = 0; i < kNumEnergyBinsFar+1; ++i ){
    //CC
    if(BMbybeams[beamToUse]->GetMCHists().at(kCCPurity_FD)->GetBinContent(i) > 0 ){
      BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPred_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPurCorr_FD)->GetBinContent(i)
                                                                                 / BMbybeams[beamToUse]->GetMCHists().at(kCCPurity_FD)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPred_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_FD)->GetBinContent(i)
                                                                                           / BMbybeams[beamToUse]->GetMCHists().at(kCCPurity_FD)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kCCRecoEnergyPred_FD)->SetBinContent(i, 0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyPred_FD)->SetBinContent(i, 0);
    }  
  
    //NC
    if(BMbybeams[beamToUse]->GetMCHists().at(kNCPurity_FD)->GetBinContent(i) > 0 ){
      BMbybeams[beamToUse]->GetMCHists().at(kNCRecoEnergyPred_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetMCHists().at(kNCRecoEnergyPurCorr_FD)->GetBinContent(i)
                                                                                 / BMbybeams[beamToUse]->GetMCHists().at(kNCPurity_FD)->GetBinContent(i));
      BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyPred_FD)->SetBinContent(i, BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyPurCorr_FD)->GetBinContent(i)
                                                                                           / BMbybeams[beamToUse]->GetMCHists().at(kNCPurity_FD)->GetBinContent(i));
    }
    else{
      BMbybeams[beamToUse]->GetMCHists().at(kNCRecoEnergyPred_FD)->SetBinContent(i, 0);
      BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyPred_FD)->SetBinContent(i, 0);
    }  
  
  }
  return;  
}

void NCExtrapolationBeamMatrix::
DoOscillations(const NC::OscProb::OscPars* pars, int beamToUse){

  
  //means a vector of arrays!
  static vector<double*> survivalProbability;
  
  //means a basic vector of doubles!
  vector<double> OscBeamMatrixCC(kNumEnergyBinsFar+1, 0.);  
  vector<double> OscBeamMatrixNC(kNumEnergyBinsFar+1, 0.);  

  static bool first = true;
  if(first){
    first = false;
    
    fNominalBM.resize(BMbybeams[beamToUse]->Get2DOscHists().size());
    fOscBM.resize(BMbybeams[beamToUse]->Get2DOscHists().size());
    survivalProbability.resize(BMbybeams[beamToUse]->Get2DOscHists().size());    
    
    for(unsigned int h = 0; h < BMbybeams[beamToUse]->Get2DOscHists().size(); ++h){
      fNominalBM[h].resize(kNumEnergyBinsFar);
      fOscBM[h].resize(kNumEnergyBinsFar);
      survivalProbability[h] = new double[kNumTrueEnergyBins];
    }
  }
    
  
  static const NC::OscProb::OscPars* prevpars = 0;
  
  // Transition probability does not depend on reco energy - so we
  // can precalculate it out here.
  // Reuse the calculations from last time round if the oscillations
  // parameters haven't changed since then.
  if(!prevpars || !pars->IsEquiv(prevpars)){
    for(int h = 0; h < int(BMbybeams[beamToUse]->Get2DOscHists().size()); ++h){
      NCType::EOscMode oscType;
      NCType::EEventType interaction;
      
      switch(h){
      case kNCSig:
      case kCCBkg:
      case kNCWrongSignBkg:
        interaction = NCType::kNC;
        oscType     = NCType::kNuMuToNuMu;
        break;
      case kCCSig:
      case kNCBkg:
      case kCCWrongSignBkg:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuMu;
        break;
      case kNCTau:
      case kCCTau:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuTau;
        break;
      case kNCBeamNue:
      case kCCBeamNue:
        interaction = NCType::kCC;
        oscType     = NCType::kNuEToNuE;
        break;
      case kNCOscNue:
      case kCCOscNue:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuE;
        break;
      default:
        MSG("NCExtrapolationBeamMatrix", Msg::kError) << "Interaction type or flavor is incorrect" << endl;
        return;
      }

      double* sph = survivalProbability[h];
      for(int j = 0; j < kNumTrueEnergyBins; ++j){
        sph[j] = pars->TransitionProbability(oscType,
                                             interaction,
                                             NCType::kBaseLineFar,
                                             (true_bins[j+1] + true_bins[j]) /2.0 );
      } // end for j
    } // end for h
  } // end if
  
  delete prevpars;
  // Sadly OscPars::Clone isn't labelled const for technical reasons.
  // Need to cast away constness so that we can call it.
  prevpars = const_cast<NC::OscProb::OscPars*>(pars)->Clone();
  
  
  //This is the original simple version
  /*
  for(int i = 0; i < kNumEnergyBinsFar; ++i ){
    for(int h = 0; h < int(BMbybeams[beamToUse]->Get2DOscHists().size()); ++h){
      fNominalBM[h][i] = 0.;
      fOscBM[h][i] = 0.;
      for(int j = 0; j < kNumTrueEnergyBins; ++j){
        //double xsecFitParam = mScaleNC.at(WhichFitParam(j+1)) ;
        double xsecFitParam = 1.0 ;
        double jbineq = (true_bins[j+1] + true_bins[j]) /2.0 ; //(j+1)*0.1 ;// make binwidth
        //if(jbineq > 120.0) jbineq = 120.0 ;
        //if(h == kNCSig ) MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "PURITAN BIN j " << j << " jbineq " << jbineq << " mScaleNC.at(WhichFitParam(jbineq)) " << mScaleNC.at(WhichFitParam(jbineq)) << endl;
        if(h == kNCSig || h == kNCWrongSignBkg || h == kCCBkg) xsecFitParam = mScaleNC.at(WhichFitParam(jbineq)) ;  
        const double tmpd = BMbybeams[beamToUse]->Get2DOscHists()[h]->GetCellContent( i+1, j+1 );
        fNominalBM[h][i] += tmpd; // *xsecFitParam;
        fOscBM[h][i] += survivalProbability[j][h]*tmpd*xsecFitParam;;
      }//end loop over true bins
    }//end loop over histograms

  
    //MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "BIN i " << i-1 
    //                                  << " fOscBM[NC][i] " << fOscBM[0][i] 
    //                                  << " fOscBM[CC][i] " << fOscBM[1][i]
    //                                  << endl ;    
  }
  */
  // And this is the optimized version
  const int nx = kNumEnergyBinsFar+2;
  for(unsigned int h = 0; h < BMbybeams[beamToUse]->Get2DOscHists().size(); ++h){
    const double* t2ra =  BMbybeams[beamToUse]->Get2DOscHists()[h]->fArray;
    vector<double>& nh = fNominalBM[h];
    vector<double>& oh = fOscBM[h];
    for(int i = 0; i < kNumEnergyBinsFar; ++i){
      double& nhi = nh[i];
      double& ohi = oh[i];
      nhi = 0;
      ohi = 0;
      const double* sph = survivalProbability[h];
      const int oset = i+1+nx;
      for(int j = 0; j < kNumTrueEnergyBins; ++j){
        double xsecFitParam = 1.0 ;
        double jbineq = (true_bins[j+1] + true_bins[j]) /2.0 ;
        //using the nd xsecweights on the Wrong sign bkg and the nominal NC should you wish to study MC only
        if(h == kNCSig || h == kNCWrongSignBkg) xsecFitParam = BMbybeams[beamToUse]->GetmScaleNC().at(WhichFitParam(jbineq)) ;  
        //if(h == kNCSig || h == kNCWrongSignBkg || h == kCCBkg) xsecFitParam = BMbybeams[beamToUse]->GetmScaleNC().at(WhichFitParam(jbineq)) ;  
        const double tmpd = t2ra[nx*j+oset];
        nhi += tmpd;
        ohi += sph[j]*tmpd*xsecFitParam;
      }//end loop over true bins 
    } // end for i
  }//end loop over histograms


  //need to truetoreco the BM prediction.
  //4 step plan 
  //1. multiply oscs into true  
  //2. rebin 2d histogram to be farBinning by farBinning
  //3. normalise new 2d histogram
  //4. true to reco
  
  vector<double> tempOscCC(kNumEnergyBinsFar+1,0) ;
  vector<double> tempOscNC(kNumEnergyBinsFar+1,0) ;

  //need to use nominal for normalisation or Systematics cancel at true to reco stage
  const NCBeam::Info beamInfoLE = NCBeam::Info(BeamType::kL010z185i, NC::RunUtil::kRunI);
  int leNomBeam = fBMbybeams.find(beamInfoLE)->second ;
  
  t2rCCRebin->Reset("ICE");
  t2rNCRebin->Reset("ICE");

  t2rCC = (TH2D*)CloneFast(BMbybeams[beamToUse]->Get2DHists()[kCCTrueRecoSig_FD]);
  t2rNC = (TH2D*)CloneFast(BMbybeams[beamToUse]->Get2DHists()[kNCTrueRecoSig_FD]);
  
  //step 1
  for(int i = 1; i < kNumTrueEnergyBins+1; ++i ){//true
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
      const double spCC = survivalProbability[kCCSig][i-1];
      const double spNC = survivalProbability[kNCSig][i-1];
      t2rCC->SetBinContent(i,j, spCC * t2rCC->GetBinContent(i,j) );
      t2rNC->SetBinContent(i,j, spNC * t2rNC->GetBinContent(i,j) );
    }
  }
  //step 2

  rebinRecoTrueToRecoRecoHist(t2rCC, t2rCCRebin); 
  rebinRecoTrueToRecoRecoHist(t2rNC, t2rNCRebin); 

  //step 3
  normaliseRecoToTrue(t2rCCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kCCTrueEnergy_FD ) );
  normaliseRecoToTrue(t2rNCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kNCTrueEnergy_FD ) );
  
  //step 4
  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//true
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
      //OscBeamMatrixCC[j-1] += ( BMbybeams[beamToUse]->GetMCHists().at(kCCTrueEnergy_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j)) ; 
      //OscBeamMatrixNC[j-1] += ( BMbybeams[beamToUse]->GetMCHists().at(kNCTrueEnergy_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j)) ; 
      
      OscBeamMatrixCC[j-1] += ( BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j)) ; 
      OscBeamMatrixNC[j-1] += ( BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j)) ; 
      
    }
  }

  BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergySig_FD)->Reset("ICE");
  BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyWrongSignBkg_FD)->Reset("ICE");
  
  BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergySig_FD)->Reset("ICE");
  BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyWrongSignBkg_FD)->Reset("ICE");

  
  for(int i = 0; i < kNumEnergyBinsFar; ++i ){
    
    // might want the MC as well at some point
    
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergySig_FD)->SetBinContent(i+1,OscBeamMatrixNC[i] );
    BMbybeams[beamToUse]->GetDataPredHists().at(kNCDataRecoEnergyWrongSignBkg_FD)->SetBinContent(i+1,fOscBM[kNCWrongSignBkg][i] );
    
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergySig_FD)->SetBinContent(i+1,OscBeamMatrixCC[i]);
    BMbybeams[beamToUse]->GetDataPredHists().at(kCCDataRecoEnergyWrongSignBkg_FD)->SetBinContent(i+1,fOscBM[kCCWrongSignBkg][i]);


    const double totalMCNCBM_fd = OscBeamMatrixNC[i] //fOscBM[kNCSig][i] //OscBeamMatrixNC[i]      
      + fOscBM[kNCWrongSignBkg][i] 
      + fOscBM[kNCBkg][i]
      + fOscBM[kNCTau][i]
      + fOscBM[kNCBeamNue][i]
      + fOscBM[kNCOscNue][i];
    
    const double totalMCCCBM_fd = OscBeamMatrixCC[i] //fOscBM[kCCSig][i] //OscBeamMatrixCC[i]    
      + fOscBM[kCCWrongSignBkg][i] 
      + fOscBM[kCCBkg][i]
      + fOscBM[kCCTau][i]
      + fOscBM[kCCBeamNue][i]
      + fOscBM[kCCOscNue][i];
    
    
    //if(fDoExtrapolation){
    BMbybeams[beamToUse]->GetFitHists()[kNCSig]->SetBinContent(i+1,totalMCNCBM_fd);
    BMbybeams[beamToUse]->GetFitHists()[kCCSig]->SetBinContent(i+1,totalMCCCBM_fd);
    //}
  }

  delete t2rNC ;
  delete t2rCC ;
  
  OscBeamMatrixNC.clear();
  OscBeamMatrixCC.clear();
}
//this will need to come AFTER the addevent loops
//don't forget to check what is scaled to what... 
void NCExtrapolationBeamMatrix::
doNDXsectionFit(int PrintLevel, NCBeam* beam_ndLE, int bemequ){

  for (Int_t i=0; i<4; i++) {
    m0ScaleNC.at(i) = 1.0;
    m0ScaleNCErr.at(i) = 0.005;
    mScaleNC.at(i) = 1.0;
    mScaleNCErr.at(i) = 0.005;
    mScaleNeg[i] = 0.0 ;
    mScalePos[i] = 0.0 ;
    corr[i] = 0.0 ;
  }
  
  // no need to do fit for all energy types LE come first
  if(beam_ndLE->GetBeamType() != BeamType::kL010z185i){ 
    for(int i = 0 ; i < 4; i++)
      BMbybeams[bemequ]->GetmScaleNC().at(i) = mScaleNC.at(i);
    return;
  }
  //will need attending to need to when using more beams copy data into systbeams
  bool doSyst = true ;
  if(bemequ == 0) doSyst = false ;

  // the way the beams vector is constructed mean this will work
  int diff = beamsNear.size()/3 ;
  MSG("NCExtrapolationBeamMatrix",Msg::kDebug) << "LE beamsNear.size() " << beamsNear.size() << " bemequ " << bemequ << " diff " << diff <<endl;
  NCBeam* beam_ndME = beamsNear.at(bemequ+diff);
  NCBeam* beam_ndHE = beamsNear.at(bemequ+(2*diff));
 
  //might need to find a better place for this
  bmFit->beamsFit.clear();
  bmFit->beamsFit.push_back(beam_ndHE);
  bmFit->beamsFit.push_back(beam_ndME);
  bmFit->beamsFit.push_back(beam_ndLE);

  // may want to pass this in as part of the registry
  // fills the histogrmas and calsluates the scales needed. 
  ScaleBeamsToSameNumberOfEvents(kNCDataHE, doSyst);

  for(UInt_t i = 0; i < (bmFit->beamsFit).size(); ++i){
    //maybe some check on nd only here    
    int data = -1 ;
    int mc = -1;
    int beamindex = -1;
    
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL010z185i){
      data = kNCDataLE; 
      mc =  kNCMCLE;  
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired LE beamtype data " << data << " mc " << mc <<endl;
      
    }
    
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL100z200i){
      data = kNCDataME; 
      mc =  kNCMCME;     
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired ME beamtype data " << data << " mc " << mc << endl;
    }
    
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL250z200i){
      data = kNCDataHE; 
      mc =  kNCMCHE; 
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired HE beamtype data " << data << " mc " << mc << endl;
    }
    beamindex = (mc-1)/2 ;   
    
    //scaling actually done here
    ND_XSectionFit.at(mc)->Scale(bmFit->mInputScale.at(mc));
    ND_XSectionFit.at(data)->Scale(bmFit->mInputScale.at(data));
    ND_XSectionFitRew.at(data)->Scale(bmFit->mInputScale.at(data));
    for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++){
      ND_XsecByParams[beamindex][ipar]->Scale(bmFit->mInputScale.at(mc));
    }
  }
      
  //do the fit
  Fit(PrintLevel);
  for(int i = 0 ; i < 4; i++)
    BMbybeams[bemequ]->SetmScaleNCparam(i,mScaleNC.at(i));
  
  MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
    << "Fitted Parameters: 0-4 " << BMbybeams[bemequ]->GetmScaleNC().at(0)  
    << " 4-8 " << BMbybeams[bemequ]->GetmScaleNC().at(1)  
    << " 8-15 " << BMbybeams[bemequ]->GetmScaleNC().at(2) 
    << " over 15 " << BMbybeams[bemequ]->GetmScaleNC().at(3) 
    << " size " << BMbybeams[bemequ]->GetmScaleNC().size() 
    << endl;
  
}
//adapted from T Raufer
void NCExtrapolationBeamMatrix::
Fit(int PrintLevel){

  Double_t arglist[10];
  Int_t    ierr;

  TMinuit* gMinuit = new TMinuit(4);

  gMinuit->SetFCN(LogLikelihoodFunc); // Likelihood fit

  arglist[0] = PrintLevel;
  gMinuit->mnexcm("SET PRInt",arglist,1,ierr);
  if ( PrintLevel<-1 ) 
    gMinuit->mnexcm("SET NOWarning",arglist,0,ierr);

  arglist[0] = 1; // added a factor 2 to the likelihood function
  gMinuit->mnexcm("SET ERRor",arglist,1,ierr);
  
  arglist[0] = 1; // minimization strategy
  gMinuit->mnexcm("SET STR",arglist,1,ierr);
  
  // Set starting values
  Double_t par[4];
  for (Int_t i=0; i<4; i++)
    par[i] = m0ScaleNC.at(i);

  // set starting error
  Double_t step[4];
  for (Int_t i=0; i<4; i++)
    step[i] = m0ScaleNCErr.at(i);

  // define parameters
  for (Int_t i=0; i<4; i++) {
    char* name = Form("scaleNC_%i",i);
    //gMinuit->DefineParameter(i,name, par[i], step[i],0., 0.); // unbounded
    gMinuit->DefineParameter(i,name, par[i], step[i],0., 5.0); // bounded
  }

  MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "Doing BOUNDED FIT 0 to 5 " << endl;
  // start restricted fit first
 
  gMinuit->FixParameter(0);
  gMinuit->FixParameter(1);
  gMinuit->FixParameter(2);

  // do a scan
 //  arglist[0] = 5; // scan parameter no 4+1
//   arglist[1] = 40;
//   arglist[2] = 0.8;
//   arglist[3] = 1.0;
//   gMinuit->mnexcm("SCAn",arglist,4,ierr);
//   TCanvas *cnew = new TCanvas("cnew","cnew");
//   TGraph *gr = (TGraph*)gMinuit->GetPlot();
//   cnew->cd();
//   gr->Draw("ap");
  
  gMinuit->Migrad();

  gMinuit->Release(2);  
  gMinuit->Migrad();

  gMinuit->Release(1);  
  gMinuit->Migrad();

  gMinuit->Release(0);  
  gMinuit->Migrad();

  // do MINOS error analysis
  gMinuit->mnmnos();

  // check convergence (still to do)
  TString convergence = gMinuit->fCstatu; 
  MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "Convergence status: "
                            << convergence << endl;

  // get result
  for (Int_t i=0; i<4; i++)
    gMinuit->GetParameter( i, mScaleNC.at(i),  mScaleNCErr.at(i) );

  //Double_t mScaleNeg[4]; // positive and negative errors 
                                          // from MINOS
  //Double_t mScalePos[4];
  //Double_t corr[4];
  // get proper asymmetric errors from MINOS fit
  for (Int_t i=0; i<4; i++)
    gMinuit->mnerrs(i,mScalePos[i], mScaleNeg[i],
                    mScaleNCErr[i],
                    corr[i] );

  // plot result  
  
  if ( PrintLevel >-2 ) {

    Double_t LL;
    Double_t pars[4];
    Int_t    npar = 4;

    for (Int_t i=0; i<4; i++)
      pars[i] = mScaleNC.at(i);
    //should return what you put in
    LogLikelihoodFunc(npar, pars, LL, pars, 0);


    MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
      << "The final likelihood is: LL=" << setprecision(10) << LL << endl;
    MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
      << "Fitted Parameters: 0-4 " << mScaleNC.at(0)  
      << " 4-8 " << mScaleNC.at(1)  
      << " 8-15 " << mScaleNC.at(2) 
      << " over 15 " << mScaleNC.at(3) << endl;

    // write some stuff to a log file
    *logFile << "Fit results:" << endl;   
    logFile->setf(ios::fixed);
    logFile->precision(5);
    for (Int_t i=0; i<4; i++) 
      *logFile << "Par " << i << ": " << mScaleNC[i] << " +"
               << mScalePos[i] << " " << mScaleNeg[i] << " +/- "
               << mScaleNCErr.at(i) << ", Corr: " << corr[i] << endl;
  
    *logFile << "Convergence Status: " << convergence << endl
             << "Final Likelihood LL=" << LL << endl;
      // << "Written histograms to " << histofile->GetName() << endl;
    logFile->unsetf(ios::fixed);
    
    /*
    //not using yet...
    // get the number of d.o.f. from chi2 fit
    Int_t dof = 0;
    if (fUseChi2Function) {
      for (UInt_t i = 0; i<fDataFraction.size(); ++i ) {
        Double_t chi2 = -1;
        Int_t ndf   -1;
        Int_t igood = 0;
        HistEvisNCd.at(i)->Chi2TestX(HistEvisNCrew.at(i),chi2,ndf,igood, 
                                     "UUCHI2");
        dof += ndf;
        MSG("NCExtrapolationRS",Msg::kDebug) << "Ndf for beam " << i << ": " 
                                             << ndf << endl; 
        if (igood != 0)
          MSG("NCExtrapolationRS",Msg::kWarning) << "Chi2 igood parameter" 
                                                 <<" is non-zer: "
                                       << igood << endl;
      }
    }
    */

    //This bit of code may be useful if one wants to write out a result tree

    // write same info to fit tree
    //    fFitRes.likelihood = LL;
    //     for (Int_t i=0; i<4; i++) {
    //       fFitRes.parNo = i;
    //       fFitRes.inputScale = mInputScale[i];
    //       fFitRes.par = (Float_t) mScaleNC[i];
    //       fFitRes.err = (Float_t) mScaleNCErr[i];
    //       fFitRes.errPos = (Float_t) mScalePos[i];
    //       fFitRes.errNeg = (Float_t) mScaleNeg[i];
    //       fFitRes.dof = dof;
    //       fFitResTree->Fill();    
    //     }
  }  

  delete gMinuit;
}
//pars is fitpars here
void NCExtrapolationBeamMatrix::
FillNDHistsForXSectionFit(NCBeam* beam_nd,
                          NCType::EEventType nccc,
                          bool doSyst){//,
                          //                      NC::SystPars systPars){

  int data = -1; 
  int mc = -1;
  int beamindex = -1;
 
  if( beam_nd->GetBeamType() == BeamType::kL010z185i){
    data = kNCDataLE; 
    mc =  kNCMCLE;     
  }
 
  if( beam_nd->GetBeamType() == BeamType::kL100z200i){
    data = kNCDataME; 
    mc =  kNCMCME;     
  }
  
  if( beam_nd->GetBeamType() == BeamType::kL250z200i){
    data = kNCDataHE; 
    mc =  kNCMCHE; 
  }
  beamindex = (mc-1)/2;
  //no cc at the moment. May not need it.

  ND_XSectionFit[data]->Reset("ICE");
  ND_XSectionFitRew[data]->Reset("ICE");
  for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++)
    ND_XsecByParams[beamindex][ipar]->Reset("ICE");
  
  ND_XSectionFit[mc]->Reset("ICE");

  for(int b  = 0; b < beam_nd->GetNumberEnergyBins(nccc); ++b){
    NCEnergyBin* bin = beam_nd->GetEnergyBin(b, nccc);

    const int dataSize       = bin->GetDataVectorSize();
    const int mcSize         = bin->GetMCSignalVectorSize();
    const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
    const int mcSize_tau     = bin->GetMCNuTauVectorSize();
    const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
    const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
    
    MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "data " << data << " data size " << dataSize  <<" mc  " << mc << " mc size" << mcSize << endl;
     

    // The variables that get filled with event information
    double trueEnergy    = 0;
    double recoShowerE   = 0;
    double recoMuonE     = 0;
    double energy        = 0;
    double weight        = 0;
    int flavor = 0;
    for(int e = 0; e < dataSize; ++e){
      bin->GetDataInformation(energy, weight, e);
      FillHistogramSmoothed(ND_XSectionFit[data], energy, weight);
      FillHistogramSmoothed(ND_XSectionFitRew[data], energy, weight);
    }
    //------------------------------
    // numu
    for(int e = 0; e < mcSize; ++e){ 
      bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(ND_XsecByParams[beamindex][WhichFitParam(trueEnergy)], recoShowerE + recoMuonE, weight);
    }

    for(int e = 0; e < mcSize_bkg; ++e){
      bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
 }

    //------------------------------
    // nutau
    for(int e = 0; e < mcSize_tau; ++e){
      bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
    }

    //------------------------------
    // beam nue
    for(int e = 0; e < mcSize_beamnue; ++e){
      bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
    }

    //------------------------------
    // oscillated nue
    for(int e = 0; e < mcSize_oscnue; ++e){
      bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
    }
  } // end for b

  
  for(int i = 1 ; i < kNumEnergyBinsFar+1 ; i++){
    double binsum = 0.0;
    for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++){
      binsum += ND_XsecByParams[beamindex][ipar]->GetBinContent(i);
    }
    
    ND_XSectionFit[mc]->SetBinContent(i,binsum);
    
    // this is a real beam and therefore has data in it. Copy to relevant dataNDHist
    if(!doSyst){  

      if(data == kNCDataLE)
        dataNDhistLE->SetBinContent(i,ND_XSectionFit[data]->GetBinContent(i));
      
      if(data == kNCDataME) 
        dataNDhistME->SetBinContent(i,ND_XSectionFit[data]->GetBinContent(i));
      
      if(data == kNCDataHE) 
        dataNDhistHE->SetBinContent(i,ND_XSectionFit[data]->GetBinContent(i));
      
    }
      
    //this is a systematically shifted beam and has no data in it. Copy from real beam
    if(doSyst){ 
      if(data == kNCDataLE){
        ND_XSectionFit[data]->SetBinContent(i,dataNDhistLE->GetBinContent(i));
        ND_XSectionFitRew[data]->SetBinContent(i,dataNDhistLE->GetBinContent(i));
          
        }      
        if(data == kNCDataME){ 
          ND_XSectionFit[data]->SetBinContent(i,dataNDhistME->GetBinContent(i));
          ND_XSectionFitRew[data]->SetBinContent(i,dataNDhistME->GetBinContent(i));
        }
        if(data == kNCDataHE){ 
          ND_XSectionFit[data]->SetBinContent(i,dataNDhistHE->GetBinContent(i));
          ND_XSectionFitRew[data]->SetBinContent(i,dataNDhistHE->GetBinContent(i));
        }
    } 
  }
}
int NCExtrapolationBeamMatrix::
WhichFitParam(double trueEnergy){
  if(trueEnergy <=4.0) return 0;
  else if(trueEnergy <= 8.0) return 1;
  else if(trueEnergy <= 15.0) return 2;
  else if(trueEnergy <= 120.0) return 3;
  // will need to put errors in at some point
  else return 4;
} 


//code ripped straight from NCExtrapolationRS.cxx
//inputs are minuit things I guess
//
void NCExtrapolationBeamMatrix::
LogLikelihoodFunc(Int_t &npar, 
                  Double_t *gin, 
                  Double_t &f, 
                  Double_t *pars, 
                  Int_t iflag)
  {  
  iflag++; // NOT USED FOR ANYTHING, JUST TO QUIET THE COMPILER
  if (0)  cout << gin << npar << endl; 

  f = 0;

//   //  cout << "Parameter 5: " << pars[4]; 
  
//   // cross section reweighting MC using neugen parameters
//   //  sFit->Reweight(ma_qe,ma_res,kno112,kno122,kno113,kno123,false);
  
// reset reweighted histograms
  for(UInt_t i = 1; i < bmFit->ND_XSectionFitRew.size(); i+=2){
    bmFit->ND_XSectionFitRew.at(i)->Reset("ICE");
  }
  
  // now fill Reweighted Histograms
  // do something faster than histogram filling
  
  for(UInt_t i = 0; i < (bmFit->beamsFit).size(); ++i){
    //maybe some check on nd only here    
    int data = -1 ;
    int mc = -1;
    int beamindex = -1;
  
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL010z185i){
      data = kNCDataLE; 
      mc =  kNCMCLE;  
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired LE beamtype data " << data << " mc " << mc << endl;
    
    }
    
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL100z200i){
      data = kNCDataME; 
      mc =  kNCMCME;     
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired ME beamtype data " << data << " mc " << mc << endl;
    }
    
    if( bmFit->beamsFit.at(i)->GetBeamType() == BeamType::kL250z200i){
      data = kNCDataHE; 
      mc =  kNCMCHE; 
      MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "acquired HE beamtype data " << data << " mc " << mc << endl;
    }
    beamindex = (mc-1)/2 ;   
    
    // loop over the 4 energy ranges in the true NC selected spectrum
    for(int ibin = 1 ; ibin < kNumEnergyBinsFar+1 ; ibin++){
      double binsum = 0.0;
      for(UInt_t ipar = 0 ; ipar <  bmFit->ND_XsecByParams.at(0).size()-1 ; ipar++){
        binsum +=  bmFit->ND_XsecByParams[beamindex][ipar]->GetBinContent(ibin) * pars[ipar];
      }
      // add in the contribution from all the backgrounds
      binsum +=  bmFit->ND_XsecByParams[beamindex][bmFit->ND_XsecByParams.at(0).size()-1]->GetBinContent(ibin);
      bmFit->ND_XSectionFitRew.at(mc)->SetBinContent(ibin,binsum);
    }
  
    
    // now calculate the LL contributions (including under and overflows
    
    // fitting with Maximum Likelihood will revisit this if I need to
    //if (!sFit->fUseChi2Function) { 
    MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "Using LL test..." << endl;
    
    // this is for fitting without the bottom two bins; again revisit if I need to...
    //if (sFit->fCutLowBins) 
    //NCStartBin +=3; // cut out bins 0-2
    // get the index right on these things...
    int startbin = 1;
    int endbin =  bmFit->ND_XSectionFitRew.at(mc)->GetNbinsX()+1 ;
    for(int ibin = startbin ; ibin < endbin; ibin++){
      if(bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin) < 0 ){
        MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "bmFit->ND_XSectionFitRew.at(mc)->GetBinContent("<<ibin<<") Gone Negative " 
                                                    << bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin) << " setting to 1.0 " << endl;
      }
      f -= bmFit->ND_XSectionFitRew.at(data)->GetBinContent(ibin)
        * log(bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin));
    }
    f +=  bmFit->ND_XSectionFitRew.at(mc)->Integral(startbin,endbin); 
  }
//    else {
// fitting with a Chi2 function
//  MAXMSG("NCExtrapolationRS",Msg::kInfo,1) << "Using Chi2 test..." << endl;
//  f += sFit->GetChiSquare(0,i); // add chi2 for NC histos  
//}
// loop over fBeams.size()
  
  f *= 2; // include factor 2 in order to keep error definitions 
  // as in a chi2 function
  // this also means, the UP value for errors is 1 (not 0.5)
  //reseting the ones that need resetting
  for(UInt_t i = 1; i < bmFit->ND_XSectionFitRew.size(); i+=2){
    bmFit->ND_XSectionFitRew.at(i)->Reset("ICE");
  }

  return;
  }
void NCExtrapolationBeamMatrix::
ScaleBeamsToSameNumberOfEvents(int beamType, bool doSyst){

  //need to scale all histograms to the same number of events. For now scale them to the HE MC stats? Its what Tobi did.
  // I guess it has the lowest stats when you get the full set.

  int other = beamType+1 ;
  if(beamType%1 == 0) other = beamType-1 ;

  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "Scaling beams to beamType: " << beamType << " " <<  other <<endl;
 
 // set pars to 1.0 to get the unscaled histograms
 //double pars[4] = {1.0,1.0,1.0,1.0};
  for(UInt_t i = 0; i < (beamsFit).size(); ++i){
    FillNDHistsForXSectionFit(beamsFit.at(i), NCType::kNC, doSyst);
  }
  
  //can't reweight them all to the same histogram. loses POT weighting for a start.
  //need to reweight each data mc pair to the same. in this case to data exposure
 
 for(UInt_t i = 0; i <  ND_XSectionFit.size() ; i+=2){

   if(ND_XSectionFit.at(i)->Integral() > 0){
     double scale =  ND_XSectionFit.at(beamType)->Integral() / ND_XSectionFit.at(i)->Integral();
     //scale *= 0.001 ;
     mInputScale.push_back(scale);
     mInputScale.push_back(scale);
   }
   else{  
     mInputScale.push_back(0);
     mInputScale.push_back(0);
   }

}
 
 MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Scaling Beam Parameters  "
                                              <<  "mInputScale 0: " << mInputScale.at(0)
                                              <<  " mInputScale 1: " << mInputScale.at(1)
                                              <<  " mInputScale 2: " << mInputScale.at(2)
                                              <<  " mInputScale 3: " << mInputScale.at(3)
                                              <<  " mInputScale 4: " << mInputScale.at(4)
                                              <<  " mInputScale 5: " << mInputScale.at(5)
                                              <<  " mInputScale 6: " << mInputScale.at(6)
                                              <<  " mInputScale 7: " << mInputScale.at(7)
                                              <<  " mInputScale 8: " << mInputScale.at(8)
                                              <<  " mInputScale 9: " << mInputScale.at(9)
                                              <<  " mInputScale 10: " << mInputScale.at(10)
                                              <<  " mInputScale 11: " << mInputScale.at(11)
                                              << endl ;
}
void NCExtrapolationBeamMatrix::
rebinRecoTrueToRecoRecoHist(TH2D *t2r, TH2D *t2rRebin){

  for(int y = 1; y < kNumEnergyBinsFar+1; ++y ){
    int x = 1;
    double trueAdd = 0.0 ;
    for(UInt_t curx = 1 ; curx <  kNumTrueEnergyBins+1 ; curx++){
      if(true_bins[curx] < kEnergyBinsFar[x]){
        trueAdd += t2r->GetBinContent(curx,y); 
      }      
      else{
        trueAdd += t2r->GetBinContent(curx,y); 
        t2rRebin->SetBinContent(x,y,trueAdd);
        trueAdd = 0.;
        x++;    
      }
    }
  }
}
void NCExtrapolationBeamMatrix::
normaliseRecoToTrue(TH2D *t2rRebin, TH1D *norm){

  for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
    for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
      double trueEnergyBinValue = norm->GetBinContent(i) ;
      if( trueEnergyBinValue > 0) t2rRebin->SetBinContent(i,j, t2rRebin->GetBinContent(i,j)/trueEnergyBinValue );
      else t2rRebin->SetBinContent(i,j,0);
      
    }
  }
}

---