NuMatrixMethod.cxx

Go to the documentation of this file.

// $Id: NuMatrixMethod.cxx,v 1.33 2009/10/01 18:56:09 bckhouse Exp $
//
// Performs the Matrix Method extrapolation.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
// Adapted from code by Chris Smith
//
// $Log: NuMatrixMethod.cxx,v $
// Revision 1.33  2009/10/01 18:56:09  bckhouse
// Centralize calculation of oscillation functions into NuUtilities: function OscillationWeight, DecayWeight, DecoherenceWeight.
//
// Revision 1.32  2008/04/14 15:39:06  evans
// Making the optimised compilier happier by removing unitialised
// variable possibilities.
//
// Revision 1.31  2008/03/06 14:51:40  evans
// Making the way the (tiny) numubar background is handled more explicit.
//
// Revision 1.30  2008/02/29 12:01:26  evans
// Putting an 'if' in so the tau code is only run for the PRL-style
// analysis. No it shouldn't crash Alex's code.
//
// Revision 1.29  2008/02/28 12:30:39  evans
// Adding an underscore to a variable name.
//
// Revision 1.28  2008/02/28 11:52:55  evans
// A fix to the tau cross section correction function: it now uses tau
// cross sections.
//
// Revision 1.27  2008/02/28 02:28:59  ahimmel
//
//
// Fixed a bug in the Transition methods.  Oops.  Problem solved.
//
// Revision 1.26  2008/02/28 01:24:31  ahimmel
//
//
// Fixing (hopefully) the bias problem for the Transition analysis.  The problem was that the NuMu efficiency was not being divided off when making the PotentialApperanceFlux.  The SetFDAppearedFidFlux and GetFDPotentialAppearanceFidFlux have been overhauled so that their actions are now clear.  They are also now nearly identical to one another (as they should be) which the only differences being the sources of the initial histogram value and whether cross section and efficency are multiplied or divided.
//
// Revision 1.25  2008/02/26 14:36:46  evans
// The PRL fitter now includes CC tau neutrinos in its FD prediction. It
// inverse-oscillates the numus, corrects their cross sections to be
// those for nutaus, corrects for nutau reconstruction and selection
// efficiencies, puts them through a nutau true->reco matrix, and adds
// them into the predictions. It's all very reminiscent of the nubar
// appearance analysis.
//
// Revision 1.24  2008/02/25 13:09:47  evans
// Pulling in the tau neutrino cross sections.
//
// Revision 1.23  2008/02/20 16:07:33  evans
// The histogram binning is now configured from the .xml files.
//
// Revision 1.22  2008/02/12 14:40:32  evans
// Correcting the decoherence function.
//
// Revision 1.21  2008/02/12 10:00:34  evans
// Adding alternative disappearance models, controled with a new argument
// in the NuMatrixFitter fitter functions (defaults to oscillations so
// existing scripts will still work).
//
// Revision 1.20  2008/02/08 16:10:19  evans
// Memory management fixes to ensure no seg faults for any analyses.
//
// Revision 1.19  2008/01/24 19:29:30  ahimmel
//
// Committing some small changes that got overlooked.  Changed NuMatrixInput so it inherits from TObject and can be stored in files.  Updated the fitters so they look for an already existing NuMatrixInput before generating one based on the file.  Fixed some kFatal messages in NuMatrixMethod.cxx.
//
// Revision 1.18  2008/01/20 01:08:19  evans
// A new setting to perform the PRL-style CC analysis.
//
// Revision 1.17  2008/01/18 07:42:15  evans
// Lots of bug fixes to get the no charge cut CC fit working.
//
// Revision 1.16  2008/01/17 01:24:26  evans
// Object can now have a NuXMLConfig passed into its constructor. If it
// does, then it does an online systematic shift to the cross section
// file it's been passed.
//
// Revision 1.15  2008/01/11 20:17:16  evans
// Changing TH1D::Clear() to TH1D::Reset().
//
// Revision 1.14  2008/01/09 01:01:19  evans
// New functions to perform the modularised no-charge-cut CC analysis in
// the way the modularised CPT analysis with fixed dm2,sn2 works, to
// allow systematics studies to be performed on the no-charge-cut CC
// analysis with the existing mechanism.
//
// Revision 1.13  2008/01/05 23:59:07  ahimmel
//
//
// - Added a function WriteInputForFitter to write a NuMatrixInput file with the results of the extrapolation.  This writes the same information that is in the WriteFileForFitter but in a different format.
//
// - Added an fextrpScheme variable so the object keeps track of what extrapolation scheme it is using.
//
// Revision 1.12  2007/12/24 11:51:02  evans
// Adding a modular fitter for the numu->numubar transition analysis. The
// interface is identical to the CPT version. This has required some
// additions to the NuMatrixOutput class (holding an extra histogram and
// a best fit transition probability), so the ClassDef has been
// incremented. But only additions have been made: no data members have
// been removed.
//
// Revision 1.11  2007/12/12 10:06:57  evans
// Allowing the cross section stage to be switched off (for a combined
// numu+numubar extrapolation in which the cross sections are included in
// the beam matrix).
//
// Revision 1.10  2007/12/05 12:31:33  evans
// Updating the numu->numubar transition analysis so that no oscillation
// parameters are required until the final modular fit stage.
//
// Revision 1.9  2007/11/30 16:40:08  evans
// Adding code to nu a numu->numubar appearance analysis with the matrix
// method.
//
// Revision 1.8  2007/11/20 12:41:18  evans
// Debugging the modular MM fitter.
//
// Revision 1.7  2007/11/19 16:27:43  evans
// Writing a function to do a nubar CPT modular MM fit for systematics
// studies.
//
// Revision 1.6  2007/11/06 18:34:38  hartnell
//
// Modified by Justin and I to output the set of files needed for an
// external fitter to fit for dm2bar, etc.
//
// Revision 1.5  2007/09/12 14:39:42  evans
// The extrapolation and fitting is now working (likelihood fit only: I
// need to rebin the histograms before a chi2 fit will work). The code
// can either do the standard CC analysis's extrapolation, or a CPT fit.
//
// Revision 1.4  2007/09/11 17:08:17  evans
// Bug fixes. Mainly memory management in CombineMMHelpers.C. Also a problem
// with units in the oscillation formulae in NuMatrixMethod. Got rid of unused
// parameters in NuMatrixFitter.
//
// Revision 1.3  2007/09/11 13:43:01  evans
// Made more extrapolation parameters configurable. Ensured the correct cross-
// sections are picked up for each neutrino type. Other minor bug fixes.
//
// Revision 1.2  2007/09/06 12:19:30  evans
// Matrix extrapolation updated to include a CPT antineutrino extrapolation.
//
// NuMatrixMethod extrapolates one ND data spectrum to the FD. NuMatrixFitter
// is a wrapper that creates as many NuMatrixMethod objects as necessary to
// extrapolate the required individual fluxes (e.g. numus and nubars
// independently). So far two CPT analysis methods are implemented: one
// performs the standard CC extrapolation on the numus and numubars separately
// with extra stages in the FD to oscillate the CC background. The other uses
// the CC numu extrapolation to calculate the CC numu background to the CC
// numubar prediction, and vice-versa.
//
// The code still requires testing.
//
// Revision 1.1  2007/09/04 14:45:58  evans
// Begun porting the matrix method extrapolation code into NtupleUtils.
//
// NuMatrixMethod is the new class that does the extrapolation.
// RunMMPrediction.C is an example of how to run it.
//
// At present, the code only runs the standard CC extrapolation on neutrinos.
//
//

#include "TArrayD.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"

#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuMatrixInput.h"
#include "NtupleUtils/NuMatrixMethod.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuXMLConfig.h"

using std::string;
using namespace NuMMExtrapolation;

ClassImp(NuMatrixMethod)

CVSID("$Id: NuMatrixMethod.cxx,v 1.33 2009/10/01 18:56:09 bckhouse Exp $");

//____________________________________________________________________72
NuMatrixMethod::NuMatrixMethod()
  : fDoXSecStep(true),
    fNDPoT(1.0),
    fFDPoT(1.0),
    fFluxPoT(1.0),
    fNDFidMass(1.0),
    fFDFidMass(1.0),
    fFDDistance(735.0*Munits::km)
{
  fxmlConfig = 0;
  
  const NuUtilities utils;
  const vector<Double_t> recoBins = utils.RecoBins(NuBinningScheme::kUnknown);
  fNRecoBins = recoBins.size()-1;
  fRecoBinEdges = new Double_t[fNRecoBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = recoBins.begin();
         itBin != recoBins.end();
         ++itBin, ++i){
      fRecoBinEdges[i] = *itBin;
    }
  }
  
  const vector<Double_t> trueBins = utils.TrueBins(NuBinningScheme::kUnknown);
  fNTrueBins = trueBins.size()-1;
  fTrueBinEdges = new Double_t[fNTrueBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = trueBins.begin();
         itBin != trueBins.end();
         ++itBin, ++i){
      fTrueBinEdges[i] = *itBin;
    }
  }

  fRecoEnergyPurCor_ND = new
    TH1D("RecoEnergyPurCor_ND",
         "ND Reconstructed Energy After Purity Correction",
         fNRecoBins,fRecoBinEdges);
  fTrueEnergyPurCor_ND = new
    TH1D("TrueEnergyPurCor_ND",
         "ND True Energy After Purity Correction",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_ND = new
    TH1D("TrueEnergyEffCor_ND",
         "ND True Energy After Purity & Efficiency Corrections",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyFlux_ND = new
    TH1D("TrueEnergyFlux_ND",
         "ND True Energy After XSec Corrections (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyMatrix_FD = new
    TH1D("TrueEnergyMatrix_FD",
         "FD True Energy After Matrix Application (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyCC_FD = new
    TH1D("TrueEnergyCC_FD",
         "FD True CC Energy After XSec Application",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_FD = new
    TH1D("TrueEnergyEffCor_FD",
         "FD True Energy After Efficiency Correction",
         fNTrueBins,fTrueBinEdges);
  //Two special case histos for FD no oscillations:
  fRecoEnergyNoOsc_FD = new
    TH1D("RecoEnergyNoOsc_FD",
         "FD CC Reco Energy Without Oscillations (no background)",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyNoOscPred_FD = new
    TH1D("RecoEnergyNoOscPred_FD",
         "FD Predicted Reco Energy Without Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoUnoscNCBackground_FD=new 
    TH1D("UnoscNCBackground_FD",
         "Reco energy unoscillated NC background prediction",
         fNRecoBins,fRecoBinEdges);
  //After the digression:
  fTrueEnergyOsc_FD = new
    TH1D("TrueEnergyOsc_FD",
         "FD True Energy With Oscillations",
         fNTrueBins,fTrueBinEdges);
  fRecoEnergyOsc_FD = new
    TH1D("RecoEnergyOsc_FD",
         "FD Reco Energy With Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyPred_FD = new
    TH1D("RecoEnergyPred_FD",
         "FD Predicted Reco Energy",
         fNRecoBins,fRecoBinEdges);
  
  fRecoEnergyMeas_ND = 0;
  fRecoEnergyMeas_FD = 0;
  fEffCorTau_FD = 0;

  fPurity_ND = 0;          
  fRecoVsTrueEnergy_ND = 0;
  fEfficiency_ND = 0;      
  fXSec_CC_Graph = 0;
  fTau_XSec_CC_Graph = 0;
  fFDVsNDMatrixRW = 0;     
  fEfficiency_FD = 0;
  fRecoVsTrueEnergy_FD = 0;
  fPurity_FD = 0;          
  fCCContamination_FD = 0;
  fNCContamination_FD = 0;
  fRecoToTrueCCContamination_FD = 0;
  fTrueToRecoCCContamination_FD = 0;
  fSuppliedTrueUnoscCCBackground_FD = 0;
  fOtherEfficiency_FD = 0;

  fAppearedTrueEffCor_FD = new
    TH1D("AppearedTrueEffCor_FD",
         "FD appearance spectrum, true energy, efficiency corrected",
         fNTrueBins,fTrueBinEdges);

  fEfficiencyTau_FD = 0;
  fRecoVsTrueEnergyTau_FD = 0;
}

//____________________________________________________________________72
NuMatrixMethod::NuMatrixMethod(const Int_t NbinsX,
                               const Double_t LowX,
                               const Double_t UppX,
                               const Double_t NDPoT,
                               const Double_t FDPoT,
                               const Double_t FluxPoT,
                               const Double_t NDFidMass,
                               const Double_t FDFidMass,
                               const string helperfilename,
                               const string xsecfilename,
                               const string NDDatafilename,
                               const string outputFileName,
                               const NuXMLConfig* xmlConfig)
  : fDoXSecStep(true),
    fFDDistance(735.0*Munits::km)
{
  fxmlConfig = xmlConfig;
  fOutName = outputFileName;
  fhelperfilename = helperfilename;
  fxsecfilename = xsecfilename;
  fndDataFilename = NDDatafilename;

  fNRecoBins = NbinsX;
  fNTrueBins = NbinsX;
  Double_t bwidth = (UppX-LowX)/Double_t(NbinsX);
  fRecoBinEdges=new Double_t[fNRecoBins+1];
  fTrueBinEdges=new Double_t[fNTrueBins+1];
  for(Int_t i=0;i<NbinsX+1;i++){
    fRecoBinEdges[i] = LowX + i*bwidth;
    fTrueBinEdges[i] = LowX + i*bwidth;
  }

  fNDPoT = NDPoT;
  fFDPoT = FDPoT;
  fFluxPoT = FluxPoT;
  fNDFidMass = NDFidMass;
  fFDFidMass = FDFidMass;

  fRecoEnergyPurCor_ND = new
    TH1D("RecoEnergyPurCor_ND",
         "ND Reconstructed Energy After Purity Correction",
         fNRecoBins,fRecoBinEdges);
  fTrueEnergyPurCor_ND = new
    TH1D("TrueEnergyPurCor_ND",
         "ND True Energy After Purity Correction",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_ND = new
    TH1D("TrueEnergyEffCor_ND",
         "ND True Energy After Purity & Efficiency Corrections",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyFlux_ND = new
    TH1D("TrueEnergyFlux_ND",
         "ND True Energy After XSec Corrections (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyMatrix_FD = new
    TH1D("TrueEnergyMatrix_FD",
         "FD True Energy After Matrix Application (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyCC_FD = new
    TH1D("TrueEnergyCC_FD",
         "FD True CC Energy After XSec Application",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_FD = new
    TH1D("TrueEnergyEffCor_FD",
         "FD True Energy After Efficiency Correction",
         fNTrueBins,fTrueBinEdges);
  //Two special case histos for FD no oscillations:
  fRecoEnergyNoOsc_FD = new
    TH1D("RecoEnergyNoOsc_FD",
         "FD CC Reco Energy Without Oscillations (no background)",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyNoOscPred_FD = new
    TH1D("RecoEnergyNoOscPred_FD",
         "FD Predicted Reco Energy Without Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoUnoscNCBackground_FD=new
    TH1D("UnoscNCBackground_FD",
         "Reco energy unoscillated NC background prediction",
         fNRecoBins,fRecoBinEdges);
  //After the digression:
  fTrueEnergyOsc_FD = new
    TH1D("TrueEnergyOsc_FD",
         "FD True Energy With Oscillations",
         fNTrueBins,fTrueBinEdges);
  fRecoEnergyOsc_FD = new
    TH1D("RecoEnergyOsc_FD",
         "FD Reco Energy With Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyPred_FD = new
    TH1D("RecoEnergyPred_FD",
         "FD Predicted Reco Energy",
         fNRecoBins,fRecoBinEdges);
  fEffCorTau_FD = 0;
  /*
  //Get cross-sections
  TFile *xsecfile = new TFile(xsecfilename.c_str(),"READ");
  TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
  XSec_CC->SetDirectory(0);
  xsecfile->Close();
  if (xsecfile){delete xsecfile; xsecfile = 0;}
  Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
  Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
  for(int i=0;i<XSec_CC->GetNbinsX();i++) {
    x[i] = XSec_CC->GetBinCenter(i+1);
    y[i] = XSec_CC->GetBinContent(i+1);
  }
  fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}
  
  //Get the helper histograms:
  TFile* helperFile = new TFile(helperfilename.c_str(),"READ");
  fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergy_ND");
  fRecoVsTrueEnergy_ND->SetDirectory(0);
  fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergy_FD");
  fRecoVsTrueEnergy_FD->SetDirectory(0);
  
  fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRW");
  if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
  
  fEfficiency_ND = (TH1D*) helperFile->Get("Efficiency_ND");
  fEfficiency_ND->SetDirectory(0);
  fEfficiency_FD = (TH1D*) helperFile->Get("Efficiency_FD");
  fEfficiency_FD->SetDirectory(0);
  fPurity_ND = (TH1D*) helperFile->Get("Purity_ND");
  fPurity_ND->SetDirectory(0);
  fPurity_FD = (TH1D*) helperFile->Get("Purity_FD");
  fPurity_FD->SetDirectory(0);
  */
  fCCContamination_FD = 0;
  fNCContamination_FD = 0;
  fRecoToTrueCCContamination_FD = 0;
  fTrueToRecoCCContamination_FD = 0;
  fSuppliedTrueUnoscCCBackground_FD = 0;
  fOtherEfficiency_FD = 0;

  fAppearedTrueEffCor_FD = new
    TH1D("AppearedTrueEffCor_FD",
         "FD appearance spectrum, true energy, efficiency corrected",
         fNTrueBins,fTrueBinEdges);

  fEfficiencyTau_FD = 0;
  fRecoVsTrueEnergyTau_FD = 0;

  /*
  helperFile->Close();
  if (helperFile) {delete helperFile; helperFile = 0;}

  //Get the data:
  TFile* ndDataFile = new TFile(NDDatafilename.c_str(), "READ");
  fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergy_ND");
  fRecoEnergyMeas_ND->SetDirectory(0);
  ndDataFile->Close();
  if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
  */
  /*
  TFile* fdDataFile = new TFile(FDDatafilename.c_str(), "READ");
  fRecoEnergyMeas_FD = (TH1D*) fdDataFile->Get("RecoEnergy_FD");
  fRecoEnergyMeas_FD->SetDirectory(0);
  fdDataFile->Close();
  if (fdDataFile){delete fdDataFile; fdDataFile = 0;}
  */
}

//____________________________________________________________________72
NuMatrixMethod::NuMatrixMethod
(const NuBinningScheme::NuBinningScheme_t binningScheme,
 const Double_t NDPoT,
 const Double_t FDPoT,
 const Double_t FluxPoT,
 const Double_t NDFidMass,
 const Double_t FDFidMass,
 const string helperfilename,
 const string xsecfilename,
 const string NDDatafilename,
 const string outputFileName,
 const NuXMLConfig* xmlConfig)
  : fDoXSecStep(true),
    fFDDistance(735.0*Munits::km)
{
  fxmlConfig = xmlConfig;
  fOutName = outputFileName;
  fhelperfilename = helperfilename;
  fxsecfilename = xsecfilename;
  fndDataFilename = NDDatafilename;
  
  const NuUtilities utils;
  const vector<Double_t> recoBins = utils.RecoBins(binningScheme);
  fNRecoBins = recoBins.size()-1;
  fRecoBinEdges = new Double_t[fNRecoBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = recoBins.begin();
         itBin != recoBins.end();
         ++itBin, ++i){
      fRecoBinEdges[i] = *itBin;
    }
  }
  
  const vector<Double_t> trueBins = utils.TrueBins(binningScheme);
  fNTrueBins = trueBins.size()-1;
  fTrueBinEdges = new Double_t[fNTrueBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = trueBins.begin();
         itBin != trueBins.end();
         ++itBin, ++i){
      fTrueBinEdges[i] = *itBin;
    }
  }
  
  fNDPoT = NDPoT;
  fFDPoT = FDPoT;
  fFluxPoT = FluxPoT;
  fNDFidMass = NDFidMass;
  fFDFidMass = FDFidMass;

  fRecoEnergyPurCor_ND = new
    TH1D("RecoEnergyPurCor_ND",
         "ND Reconstructed Energy After Purity Correction",
         fNRecoBins,fRecoBinEdges);
  fTrueEnergyPurCor_ND = new
    TH1D("TrueEnergyPurCor_ND",
         "ND True Energy After Purity Correction",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_ND = new
    TH1D("TrueEnergyEffCor_ND",
         "ND True Energy After Purity & Efficiency Corrections",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyFlux_ND = new
    TH1D("TrueEnergyFlux_ND",
         "ND True Energy After XSec Corrections (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyMatrix_FD = new
    TH1D("TrueEnergyMatrix_FD",
         "FD True Energy After Matrix Application (flux)",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyCC_FD = new
    TH1D("TrueEnergyCC_FD",
         "FD True CC Energy After XSec Application",
         fNTrueBins,fTrueBinEdges);
  fTrueEnergyEffCor_FD = new
    TH1D("TrueEnergyEffCor_FD",
         "FD True Energy After Efficiency Correction",
         fNTrueBins,fTrueBinEdges);
  //Two special case histos for FD no oscillations:
  fRecoEnergyNoOsc_FD = new
    TH1D("RecoEnergyNoOsc_FD",
         "FD CC Reco Energy Without Oscillations (no background)",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyNoOscPred_FD = new
    TH1D("RecoEnergyNoOscPred_FD",
         "FD Predicted Reco Energy Without Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoUnoscNCBackground_FD=new
    TH1D("UnoscNCBackground_FD",
         "Reco energy unoscillated NC background prediction",
         fNRecoBins,fRecoBinEdges);
  //After the digression:
  fTrueEnergyOsc_FD = new
    TH1D("TrueEnergyOsc_FD",
         "FD True Energy With Oscillations",
         fNTrueBins,fTrueBinEdges);
  fRecoEnergyOsc_FD = new
    TH1D("RecoEnergyOsc_FD",
         "FD Reco Energy With Oscillations",
         fNRecoBins,fRecoBinEdges);
  fRecoEnergyPred_FD = new
    TH1D("RecoEnergyPred_FD",
         "FD Predicted Reco Energy",
         fNRecoBins,fRecoBinEdges);
  fEffCorTau_FD = 0;

  fCCContamination_FD = 0;
  fNCContamination_FD = 0;
  fRecoToTrueCCContamination_FD = 0;
  fTrueToRecoCCContamination_FD = 0;
  fSuppliedTrueUnoscCCBackground_FD = 0;
  fOtherEfficiency_FD = 0;

  fAppearedTrueEffCor_FD = new
    TH1D("AppearedTrueEffCor_FD",
         "FD appearance spectrum, true energy, efficiency corrected",
         fNTrueBins,fTrueBinEdges);

  fEfficiencyTau_FD = 0;
  fRecoVsTrueEnergyTau_FD = 0;
}

//____________________________________________________________________72
NuMatrixMethod
::NuMatrixMethod(const NuMatrixInput& mmInput,
                 const NuMMExtrapolation::NuMMScheme_t extrpScheme)
  : fDoXSecStep(true),
    fFDDistance(735.0*Munits::km)
{
  fxmlConfig = 0;//Shouldn't need this for FD only steps.
  if (NuMMExtrapolation::kModularNuMuCPTFit == extrpScheme){
    //M numu(bar)
    fRecoVsTrueEnergy_FD = new
      TH2D(*mmInput.NuMuTrueToRecoFD());
    //M-squigle numu(bar) for WS
    fTrueToRecoCCContamination_FD = new
      TH2D(*mmInput.NuMuTrueToRecoCCBkgFD());
    //T
    fTrueEnergyEffCor_FD = new
      TH1D(*mmInput.NuMuTrueEnFD());
    //K
    fSuppliedTrueUnoscCCBackground_FD = new
      TH1D(*mmInput.NuMuTrueEnCCBkgFD());
    //Z
    fRecoUnoscNCBackground_FD = new
      TH1D(*mmInput.NuMuRecoEnNCBkgFD());
    //numubar->numu appearance spectrum
    fAppearedTrueEffCor_FD = new
      TH1D(*mmInput.NuMuTrueEnPotentialAppearanceEffCorFD());
  }
  else if (NuMMExtrapolation::kModularNuMuBarCPTFit == extrpScheme){
    //M numu(bar)
    fRecoVsTrueEnergy_FD = new
      TH2D(*mmInput.NuMuBarTrueToRecoFD());
    //M-squigle numu(bar) for WS
    fTrueToRecoCCContamination_FD = new
      TH2D(*mmInput.NuMuBarTrueToRecoCCBkgFD());
    //T
    fTrueEnergyEffCor_FD = new
      TH1D(*mmInput.NuMuBarTrueEnFD());
    //K
    fSuppliedTrueUnoscCCBackground_FD = new
      TH1D(*mmInput.NuMuBarTrueEnCCBkgFD());
    //Z
    fRecoUnoscNCBackground_FD = new
      TH1D(*mmInput.NuMuBarRecoEnNCBkgFD());
    //numu->numubar appearance spectrum
    fAppearedTrueEffCor_FD = new
      TH1D(*mmInput.NuMuBarTrueEnPotentialAppearanceEffCorFD());
  }
  else if (NuMMExtrapolation::kModularNuMuBarTransitionFit == extrpScheme){
    //M numu(bar)
    fRecoVsTrueEnergy_FD = new
      TH2D(*mmInput.NuMuBarTrueToRecoFD());
    //M-squigle numu(bar) for WS
    fTrueToRecoCCContamination_FD = new
      TH2D(*mmInput.NuMuBarTrueToRecoCCBkgFD());
    //T
    fTrueEnergyEffCor_FD = new
      TH1D(*mmInput.NuMuBarTrueEnFD());
    //K
    fSuppliedTrueUnoscCCBackground_FD = new
      TH1D(*mmInput.NuMuBarTrueEnCCBkgFD());
    //Z
    fRecoUnoscNCBackground_FD = new
      TH1D(*mmInput.NuMuBarRecoEnNCBkgFD());
    //numu->numubar appearance spectrum
    fAppearedTrueEffCor_FD = new
      TH1D(*mmInput.NuMuBarTrueEnPotentialAppearanceEffCorFD());
  }
  else if (NuMMExtrapolation::kModularNoChargeCutFit == extrpScheme){
    //M numu(bar)
    fRecoVsTrueEnergy_FD = new
      TH2D(*mmInput.NoChargeCutTrueToRecoFD());
//     //M-squigle numu(bar) for WS
//     fTrueToRecoCCContamination_FD = new
//       TH2D(*mmInput.NoChargeCutTrueToRecoCCBkgFD());
    //T
    fTrueEnergyEffCor_FD = new
      TH1D(*mmInput.NoChargeCutTrueEnFD());
//     //K
//     fSuppliedTrueUnoscCCBackground_FD = new
//       TH1D(*mmInput.NoChargeCutTrueEnCCBkgFD());
    //Z
    fRecoUnoscNCBackground_FD = new
      TH1D(*mmInput.NoChargeCutRecoEnNCBkgFD());
    //numu->numubar appearance spectrum
//     fAppearedTrueEffCor_FD = new
//       TH1D(*mmInput.NuMuBarTrueEnPotentialAppearanceEffCorFD());
    fAppearedTrueEffCor_FD = 0;
  }
  else if (NuMMExtrapolation::kModularPRLCCFit == extrpScheme){
    //M numu(bar)
    fRecoVsTrueEnergy_FD = new
      TH2D(*mmInput.NuMuTrueToRecoFD());
    //T
    fTrueEnergyEffCor_FD = new
      TH1D(*mmInput.NuMuTrueEnFD());
    //Z
    fRecoUnoscNCBackground_FD = new
      TH1D(*mmInput.NuMuRecoEnNCBkgFD());
    fAppearedTrueEffCor_FD = 0;
    fEffCorTau_FD = new
      TH1D(*mmInput.NuMuEffCorTaus_FD());
    fRecoVsTrueEnergyTau_FD = new
      TH2D(*mmInput.NuMuRecoVsTrueTaus_FD());
  }
  else{
    try {
      MSG("NuMatrixMethod.cxx", Msg::kFatal) 
        << "Not a valid extrapolation scheme for a modular "
        << "UK matrix method fit"
        << endl;
    }
    catch(MSGException) {
      MSG("NuMatrixMethod.cxx", Msg::kError)
        << "Not a valid extrapolation scheme for a modular "
        << "UK matrix method fit"
        << endl;
    }
    
    
  }
  //Create the other histograms we need, using existing ones to get binning
  fTrueEnergyOsc_FD = new TH1D(*fTrueEnergyEffCor_FD);
  fTrueEnergyOsc_FD->Reset();
  
  fRecoEnergyOsc_FD = new TH1D(*fRecoUnoscNCBackground_FD);
  fRecoEnergyOsc_FD->Reset();
  
  fRecoEnergyPred_FD = new TH1D(*fRecoUnoscNCBackground_FD);
  fRecoEnergyPred_FD->Reset();

  fNRecoBins = fRecoEnergyPred_FD->GetNbinsX();
  const TArrayD* recoBinEdges = fRecoEnergyPred_FD->GetXaxis()->GetXbins();
  fRecoBinEdges = new Double_t[recoBinEdges->GetSize()];
  for (Int_t i=0; i<recoBinEdges->GetSize(); ++i){
    fRecoBinEdges[i] = recoBinEdges->At(i);
  }

  fNTrueBins = fTrueEnergyEffCor_FD->GetNbinsX();
  const TArrayD* trueBinEdges = fTrueEnergyEffCor_FD->GetXaxis()->GetXbins();
  fTrueBinEdges = new Double_t[trueBinEdges->GetSize()];
  for (Int_t i=0; i<trueBinEdges->GetSize(); ++i){
    fTrueBinEdges[i] = trueBinEdges->At(i);
  }

  fOutName = "";
  
  fextrpScheme = extrpScheme;

  fEfficiencyTau_FD = 0;
  //  fRecoVsTrueEnergyTau_FD = 0;

  //Zero all the histograms we don't need.
  fRecoEnergyMeas_ND = 0;  
  fRecoEnergyMeas_FD = 0;  
  fPurity_ND = 0;          
  fRecoVsTrueEnergy_ND = 0;
  fEfficiency_ND = 0;      
  fXSec_CC_Graph = 0;
  fTau_XSec_CC_Graph = 0;
  fFDVsNDMatrixRW = 0;     
  fEfficiency_FD = 0;
  fPurity_FD = 0;
  fCCContamination_FD = 0;
  fNCContamination_FD = 0;
  fRecoToTrueCCContamination_FD = 0;
  fOtherEfficiency_FD = 0;
  fRecoEnergyPurCor_ND = 0;
  fTrueEnergyPurCor_ND = 0;
  fTrueEnergyEffCor_ND = 0;
  fTrueEnergyFlux_ND = 0;
  fTrueEnergyMatrix_FD = 0;
  fTrueEnergyCC_FD = 0;
  fRecoEnergyNoOsc_FD = 0;
  fRecoEnergyNoOscPred_FD = 0;
}

//____________________________________________________________________72
NuMatrixMethod::~NuMatrixMethod()
{
  if (""!=fOutName){
    TFile *sfile = new TFile(fOutName.c_str(),"UPDATE");
    sfile->cd();
    fRecoEnergyPurCor_ND->Write();
    fTrueEnergyPurCor_ND->Write();
    fTrueEnergyEffCor_ND->Write();
    fTrueEnergyFlux_ND->Write();
    fTrueEnergyMatrix_FD->Write();
    fTrueEnergyCC_FD->Write();
    fTrueEnergyEffCor_FD->Write();
    fTrueEnergyOsc_FD->Write();
    fRecoEnergyOsc_FD->Write();
    fRecoEnergyPred_FD->Write();  
    fRecoEnergyNoOsc_FD->Write();
    fRecoEnergyNoOscPred_FD->Write();
    
    fRecoEnergyMeas_ND->Write("NDData");
    // fRecoEnergyMeas_FD->Write("FDData");
    
    sfile->Close();
  }

  if (fRecoEnergyPurCor_ND){
    delete fRecoEnergyPurCor_ND;
    fRecoEnergyPurCor_ND = 0;
  }
  if (fTrueEnergyPurCor_ND){
    delete fTrueEnergyPurCor_ND;
    fTrueEnergyPurCor_ND = 0;
  }
  if (fTrueEnergyEffCor_ND){
    delete fTrueEnergyEffCor_ND;
    fTrueEnergyEffCor_ND = 0;
  }
  if (fTrueEnergyFlux_ND){
    delete fTrueEnergyFlux_ND;
    fTrueEnergyFlux_ND = 0;
  }
  if (fTrueEnergyMatrix_FD){
    delete fTrueEnergyMatrix_FD;
    fTrueEnergyMatrix_FD = 0;
  }
  if (fTrueEnergyCC_FD){
    delete fTrueEnergyCC_FD;
    fTrueEnergyCC_FD = 0;
  }
  if (fTrueEnergyEffCor_FD){
    delete fTrueEnergyEffCor_FD;
    fTrueEnergyEffCor_FD = 0;
  }
  if (fRecoEnergyNoOsc_FD){
    delete fRecoEnergyNoOsc_FD;
    fRecoEnergyNoOsc_FD = 0;
  }
  if (fRecoEnergyNoOscPred_FD){
    delete fRecoEnergyNoOscPred_FD;
    fRecoEnergyNoOscPred_FD = 0;
  }
  if (fEffCorTau_FD){
    delete fEffCorTau_FD;
    fEffCorTau_FD = 0;
  }
}

//____________________________________________________________________72

void NuMatrixMethod::WriteFilesForFitter(std::string sPrefix)
{
  
  //helper files
  //M numu(bar)
  fRecoVsTrueEnergy_FD->
    Write((sPrefix+"TrueToRecoFD").c_str());
  //M-squigle numu(bar) for WS
  if (fTrueToRecoCCContamination_FD){
    fTrueToRecoCCContamination_FD->
      Write((sPrefix+"TrueToRecoCCBkgFD").c_str());
  }
  
  //data or mock data
  //T
  fTrueEnergyEffCor_FD->
  Write((sPrefix+"TrueEnFD").c_str());
  //K
  if (fSuppliedTrueUnoscCCBackground_FD){
    fSuppliedTrueUnoscCCBackground_FD->
      Write((sPrefix+"TrueEnCCBkgFD").c_str());
  }
  //Z
  if (fRecoVsTrueEnergy_FD){
    fRecoUnoscNCBackground_FD->
      Write((sPrefix+"RecoEnNCBkgFD").c_str());
  }

  //Appeared spectrum for numu->nubar analysis, efficiency corrected
  if (fAppearedTrueEffCor_FD){
    fAppearedTrueEffCor_FD->
      Write((sPrefix+"TrueEnPotentialAppearanceEffCor").c_str());
  }

  //Potential taus (efficiency corrected)
  if (fEffCorTau_FD){
    fEffCorTau_FD->Write((sPrefix+"EffCorTau_FD").c_str());
  }
  //Tau reco v. true energy matrix
  if (fRecoVsTrueEnergyTau_FD){
    fRecoVsTrueEnergyTau_FD->Write((sPrefix+"RecoVsTrueEnergyTau_FD").c_str());
  }
}

//____________________________________________________________________72

void NuMatrixMethod::WriteInputForFitter(NuMatrixInput *mmInput)
{
  bool charge = false;
  if (NuMMExtrapolation::kCCNuMuNoOscBackground == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuOscBackgroundUncombined == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuOscBackgroundCombined == fextrpScheme   ||
      NuMMExtrapolation::kModularNuMuCPTFit == fextrpScheme
      )
    charge = true;
  else if (NuMMExtrapolation::kCCNuMuBarOscBackgroundUncombined == fextrpScheme ||
           NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined == fextrpScheme ||
           NuMMExtrapolation::kModularNuMuBarCPTFit == fextrpScheme ||
           NuMMExtrapolation::kModularNuMuBarTransitionFit == fextrpScheme
           )
    charge = false;
  else{
    try {
      MSG("NuMatrixMethod.cxx", Msg::kFatal) 
        << "Not a valid extrapolation scheme for a modular "
        << "UK matrix method fit"
        << endl;
    }
    catch(MSGException) {
      MSG("NuMatrixMethod.cxx", Msg::kError)
        << "Not a valid extrapolation scheme for a modular "
        << "UK matrix method fit"
        << endl;
    }
  }
  
  mmInput->M(fRecoVsTrueEnergy_FD, charge);
  mmInput->Mtilde(fTrueToRecoCCContamination_FD, charge);
  mmInput->T(fTrueEnergyEffCor_FD, charge);
  mmInput->K(fSuppliedTrueUnoscCCBackground_FD, charge);
  mmInput->Z(fRecoUnoscNCBackground_FD, charge);
  mmInput->Appear(fAppearedTrueEffCor_FD, charge);
}

//____________________________________________________________________72
void NuMatrixMethod
::SetExtrapolationScheme(const NuMMScheme_t scheme)
{
  fextrpScheme = scheme;
  switch (scheme){
  case kCCNuMuNoOscBackground:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergy_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergy_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRW");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("Efficiency_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("Efficiency_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("Purity_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("Purity_FD");
      fPurity_FD->SetDirectory(0);
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numu)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (tau)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergy_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kCCNuMuOscBackgroundUncombined:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergy_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergy_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRW");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("Efficiency_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("Efficiency_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("Purity_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("Purity_FD");
      fPurity_FD->SetDirectory(0);
      
      fCCContamination_FD = (TH1D*) helperFile->Get("CCContamination_FD");
      fCCContamination_FD->SetDirectory(0);
      fNCContamination_FD = (TH1D*) helperFile->Get("NCContamination_FD");
      fNCContamination_FD->SetDirectory(0);
      fRecoToTrueCCContamination_FD = (TH2D*)
        helperFile->Get("CCContaminationRecoVsTrue_FD");
      fRecoToTrueCCContamination_FD->SetDirectory(0);
      fTrueToRecoCCContamination_FD = new TH2D(*fRecoToTrueCCContamination_FD);
      fTrueToRecoCCContamination_FD->SetDirectory(0);
      fSuppliedTrueUnoscCCBackground_FD = 0;
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numu)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (tau)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergy_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kCCNuMuBarOscBackgroundUncombined:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergyPQ_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergyPQ_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRWPQ");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("EfficiencyPQ_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("EfficiencyPQ_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("PurityPQ_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("PurityPQ_FD");
      fPurity_FD->SetDirectory(0);
      
      fCCContamination_FD = (TH1D*) helperFile->Get("CCContaminationPQ_FD");
      fCCContamination_FD->SetDirectory(0);
      fNCContamination_FD = (TH1D*) helperFile->Get("NCContaminationPQ_FD");
      fNCContamination_FD->SetDirectory(0);
      fRecoToTrueCCContamination_FD = (TH2D*)
        helperFile->Get("CCContaminationRecoVsTruePQ_FD");
      fRecoToTrueCCContamination_FD->SetDirectory(0);
      fTrueToRecoCCContamination_FD = new TH2D(*fRecoToTrueCCContamination_FD);
      fTrueToRecoCCContamination_FD->SetDirectory(0);
      fSuppliedTrueUnoscCCBackground_FD = 0;
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numubar)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (taubar)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutaubar_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergyPQ_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kCCNuMuOscBackgroundCombined:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergy_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergy_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRW");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("Efficiency_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("Efficiency_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("Purity_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("Purity_FD");
      fPurity_FD->SetDirectory(0);
      
      fCCContamination_FD = (TH1D*) helperFile->Get("CCContamination_FD");
      fCCContamination_FD->SetDirectory(0);
      fNCContamination_FD = (TH1D*) helperFile->Get("NCContamination_FD");
      fNCContamination_FD->SetDirectory(0);
      fRecoToTrueCCContamination_FD = (TH2D*)
        helperFile->Get("CCContaminationRecoVsTrue_FD");
      fRecoToTrueCCContamination_FD->SetDirectory(0);
      fTrueToRecoCCContamination_FD = new TH2D(*fRecoToTrueCCContamination_FD);
      fTrueToRecoCCContamination_FD->SetDirectory(0);
      fSuppliedTrueUnoscCCBackground_FD = 0;
      fOtherEfficiency_FD = (TH1D*)
        helperFile->Get("OtherEfficiency_FD");
      fOtherEfficiency_FD->SetDirectory(0);
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numu)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (tau)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergy_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kCCNuMuBarOscBackgroundCombined:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergyPQ_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergyPQ_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRWPQ");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("EfficiencyPQ_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("EfficiencyPQ_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("PurityPQ_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("PurityPQ_FD");
      fPurity_FD->SetDirectory(0);
      
      fCCContamination_FD = (TH1D*) helperFile->Get("CCContaminationPQ_FD");
      fCCContamination_FD->SetDirectory(0);
      fNCContamination_FD = (TH1D*) helperFile->Get("NCContaminationPQ_FD");
      fNCContamination_FD->SetDirectory(0);
      fRecoToTrueCCContamination_FD = (TH2D*)
        helperFile->Get("CCContaminationRecoVsTruePQ_FD");
      fRecoToTrueCCContamination_FD->SetDirectory(0);
      fTrueToRecoCCContamination_FD = new TH2D(*fRecoToTrueCCContamination_FD);
      fTrueToRecoCCContamination_FD->SetDirectory(0);
      fSuppliedTrueUnoscCCBackground_FD = 0;
      fOtherEfficiency_FD = (TH1D*)
        helperFile->Get("OtherEfficiencyPQ_FD");
      fOtherEfficiency_FD->SetDirectory(0);
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numubar)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (taubar)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutaubar_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergyPQ_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kCCNoChargeCut:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergyAll_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergyAll_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixXSecRWAllCC");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = false;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("EfficiencyAll_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("EfficiencyAll_FD");
      fEfficiency_FD->SetDirectory(0);
      fPurity_ND = (TH1D*) helperFile->Get("PurityAll_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("PurityAll_FD");
      fPurity_FD->SetDirectory(0);
      
      fNCContamination_FD = new TH1D(*fPurity_FD);
      fNCContamination_FD->Reset();
      for(int i=1;i<=fNCContamination_FD->GetNbinsX();i++) {
        fNCContamination_FD->SetBinContent(i,1.0);
      }
      fNCContamination_FD->Add(fPurity_FD,-1.0);
      fNCContamination_FD->SetDirectory(0);
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numu; unused, I belive)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (tau; again: unused at present)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergyAll_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  case kPRLCCAnalysis:
    {
      //Get the helper histograms:
      TFile* helperFile = new TFile(fhelperfilename.c_str(),"READ");
      fRecoVsTrueEnergy_ND = (TH2D*) helperFile->Get("RecoVsTrueEnergy_ND");
      fRecoVsTrueEnergy_ND->SetDirectory(0);
      fRecoVsTrueEnergy_FD = (TH2D*) helperFile->Get("RecoVsTrueEnergy_FD");
      fRecoVsTrueEnergy_FD->SetDirectory(0);
      
      fFDVsNDMatrixRW = (TH2D*) helperFile->Get("FDVsNDMatrixRW");
      if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
      fDoXSecStep = true;
      
      fEfficiency_ND = (TH1D*) helperFile->Get("Efficiency_ND");
      fEfficiency_ND->SetDirectory(0);
      fEfficiency_FD = (TH1D*) helperFile->Get("Efficiency_FD");
      fEfficiency_FD->SetDirectory(0);
      //The purity histograms in the helper file remove both the ND
      //background and the wrong sign CCs: the PRL analysis doesn't
      //remove the latter.
      
       fPurity_FD = (TH1D*) helperFile->Get("NCContamination_FD");
       fPurity_FD->SetDirectory(0);
       fPurity_ND = (TH1D*) helperFile->Get("NCContamination_ND");
       fPurity_ND->SetDirectory(0);
       for (int i=1; i<=fPurity_ND->GetNbinsX(); ++i){
        fPurity_ND->SetBinContent(i, 1.0 -
                                  fPurity_ND->GetBinContent(i));
        fPurity_FD->SetBinContent(i, 1.0 -
                                  fPurity_FD->GetBinContent(i));
       }
    
       fNCContamination_FD = (TH1D*) helperFile->Get("NCContamination_FD");
       fNCContamination_FD->SetDirectory(0);
      
       /*
      fPurity_ND = (TH1D*) helperFile->Get("Purity_ND");
      fPurity_ND->SetDirectory(0);
      fPurity_FD = (TH1D*) helperFile->Get("Purity_FD");
      fPurity_FD->SetDirectory(0);
       */
      //Tau helpers
      fEfficiencyTau_FD = (TH1D*) helperFile->Get("EfficiencyTau_FD");
      fEfficiencyTau_FD->SetDirectory(0);
      fRecoVsTrueEnergyTau_FD = (TH2D*)
        helperFile->Get("RecoVsTrueEnergyTau_FD");
      fRecoVsTrueEnergyTau_FD->SetDirectory(0);
      
      helperFile->Close();
      if (helperFile) {delete helperFile; helperFile = 0;}
      
      //Get cross-sections (numu)
      TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
      XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
      Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
      for(int i=0;i<XSec_CC->GetNbinsX();i++) {
        x[i] = XSec_CC->GetBinCenter(i+1);
        y[i] = XSec_CC->GetBinContent(i+1);
      }
      fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}
      
      //Get cross-sections (tau)
      xsecfile = new TFile(fxsecfilename.c_str(),"READ");
      TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
      tau_XSec_CC->SetDirectory(0);
      xsecfile->Close();
      if (xsecfile){delete xsecfile; xsecfile = 0;}
      x = new Float_t[tau_XSec_CC->GetNbinsX()];
      y = new Float_t[tau_XSec_CC->GetNbinsX()];
      for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
        x[i] = tau_XSec_CC->GetBinCenter(i+1);
        y[i] = tau_XSec_CC->GetBinContent(i+1);
      }
      fTau_XSec_CC_Graph = new TGraph(fNTrueBins,x,y);
      if (x) {delete[] x; x = 0;}
      if (y) {delete[] y; y = 0;}

      //Get the data:
      TFile* ndDataFile = new TFile(fndDataFilename.c_str(), "READ");
      fRecoEnergyMeas_ND = (TH1D*) ndDataFile->Get("RecoEnergy_ND");
      fRecoEnergyMeas_ND->SetDirectory(0);
      ndDataFile->Close();
      if (ndDataFile){delete ndDataFile; ndDataFile = 0;}
    }
    break;
  default:
    {
    }
    break;
  }
      
  if (fRecoToTrueCCContamination_FD){
    //Normalise the FD CC contamination reco to true CC matrix
    //First loop over true axis:
    for (Int_t j=1; j<=fRecoToTrueCCContamination_FD->GetNbinsY(); ++j){
      //Find out the total number of reco events in this row:
      Double_t recoEvents = 0.0;
      for (Int_t i=1; i<=fRecoToTrueCCContamination_FD->GetNbinsX(); ++i){
        recoEvents += fRecoToTrueCCContamination_FD->GetBinContent(i,j);
      }
      //Loop again, this time normalise.
      for (Int_t i=1; i<=fRecoToTrueCCContamination_FD->GetNbinsX(); ++i){
        if (recoEvents>0){
          Double_t oldContent =
            fRecoToTrueCCContamination_FD->GetBinContent(i,j);
          Double_t oldError =
            fRecoToTrueCCContamination_FD->GetBinError(i,j);
          fRecoToTrueCCContamination_FD->SetBinContent
            (i,j,oldContent/recoEvents);
          fRecoToTrueCCContamination_FD->SetBinError
            (i,j,oldError/recoEvents);
        }
        else{
          fRecoToTrueCCContamination_FD->SetBinContent(i,j,0);
          
          fRecoToTrueCCContamination_FD->SetBinError(i,j,0);
        }
      }
    }
  }
  
  if (fTrueToRecoCCContamination_FD){ 
    //Normalise the FD CC contamination true to reco CC matrix
    //First loop over reco axis:
    for (Int_t i=1; i<=fTrueToRecoCCContamination_FD->GetNbinsX(); ++i){
      //Find out the total number of true events in this column:
      Double_t trueEvents = 0.0;
      for (Int_t j=1; j<=fTrueToRecoCCContamination_FD->GetNbinsY(); ++j){
        trueEvents += fTrueToRecoCCContamination_FD->GetBinContent(i,j);
      }
      //Loop again, this time normalise.
      for (Int_t j=1; j<=fTrueToRecoCCContamination_FD->GetNbinsY(); ++j){
        if (trueEvents>0){
          Double_t oldContent =
            fTrueToRecoCCContamination_FD->GetBinContent(i,j);
          Double_t oldError =
            fTrueToRecoCCContamination_FD->GetBinError(i,j);
          fTrueToRecoCCContamination_FD->SetBinContent
            (i,j,oldContent/trueEvents);
          fTrueToRecoCCContamination_FD->SetBinError
            (i,j,oldError/trueEvents);
        }
        else{
          fTrueToRecoCCContamination_FD->SetBinContent(i,j,0);
          
          fTrueToRecoCCContamination_FD->SetBinError(i,j,0);
        }
      }
    }
  }
//   if (fxmlConfig){
//     this->ShiftXSecGraph();
//   }
}

//____________________________________________________________________72
void NuMatrixMethod::ShiftXSecGraph()
{
  NuSystematic nuSyst(*fxmlConfig);
  NuParticle::NuParticleType_t particle = NuParticle::kUndefined;
  if (NuMMExtrapolation::kCCNuMuNoOscBackground == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuOscBackgroundUncombined == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuOscBackgroundCombined == fextrpScheme ||
      NuMMExtrapolation::kModularNuMuCPTFit == fextrpScheme ||
      NuMMExtrapolation::kModularPRLCCFit == fextrpScheme){
    particle = NuParticle::kNuMu;
  }
  else if (NuMMExtrapolation::kCCNuMuBarOscBackgroundUncombined == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined == fextrpScheme ||
      NuMMExtrapolation::kModularNuMuBarCPTFit == fextrpScheme ||
      NuMMExtrapolation::kModularNuMuBarTransitionFit == fextrpScheme){
    particle = NuParticle::kNuMuBar;
  }
  else if (NuMMExtrapolation::kCCNoChargeCut == fextrpScheme ||
           NuMMExtrapolation::kModularNoChargeCutFit == fextrpScheme){
    MSG("NuMatrixMethod.cxx",Msg::kInfo)
      << "Cross section file isn't used by NuMatrixMethod for this fitting mode"
      << endl;
    particle = NuParticle::kUndefined;
    return;
  }
  else {
    MSG("NuMatrixMethod.cxx",Msg::kWarning)
      << "I don't know what I'm extrapolaing!"
      << endl;
    particle = NuParticle::kUndefined;
  }
  if (NuParticle::kUndefined != particle){
    //Loop over graph
    for (Int_t point = 0; point < fXSec_CC_Graph->GetN(); ++point){
      //Get energy & nominal cross section
      Double_t energy = 0;
      Double_t xSec = 0;
      fXSec_CC_Graph->GetPoint(point,energy,xSec);
      xSec *= nuSyst.XSecShiftScale(energy,particle);
      fXSec_CC_Graph->SetPoint(point,energy,xSec);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::PredictFDFluxUnosc()
{

  //just zero everything
  this->ResetUnoscHistograms();

  //subtract background
  this->PurityCorrectND();

  //convert ND reco spectrum into true spectrum
  this->RecoToTrueND();

  //correct for eff. in true energy
  //numerator is true CC events reco'd in the fid vol that pass all cuts
  //denominator is true CC events with true vtx in fid vol
  //hence efficiency can be >1 if more leak-in than leak-out
  this->EfficiencyCorrectND();

  //convert to flux (divide by NDmass, NDpot, xsec (xsec only if not
  //doing combined extrapolation))
  this->XSecCorrectND();

  //APPLY THE BEAM MATRIX
  this->ExtrapolateNDToFD();

  //covert from flux to events (times by FDmass, FDpot and xsec). Only
  //if fDoXSecStep == true are the cross sections included: for a
  //combined numu+numubar extrapolation, the cross sections are
  //included in the beam matrix.
  this->XSecCorrectFD();

  //correct for eff. in true energy
  //numerator is true CC events reco'd in the fid vol that pass all cuts
  //denominator is true CC events with true vtx in fid vol
  //hence efficiency can be >1 if more leak-in than leak-out
  //The "T" spectrum to pass out to external fitter is now created
  this->EfficiencyCorrectFD();

  //Now create the Z spectrum (NC events)
  //apply true to reco matrix to get pure reco CC spectrum
  this->UnoscTrueToRecoFD();

  //divide by overall purity to get the actual unoscillated FD spectrum
  this->UnoscPurityCorrectFD();

  //take unoscillated FD spectrum and times by NC contamination
  //to get the NC spectrum
  //The "Z" spectrum to pass out to external fitter is now created
  this->CalculateNCContamination();
  
  //Calculate the true energy tau spectrum as if every numu has turned
  //into a tau (only do for PRL-style analysis so far)
  if (NuMMExtrapolation::kPRLCCAnalysis == fextrpScheme){
    this->PotentialTausTrueEnergy();
  }
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDSpectrumBackgroundNoOsc(const Double_t dm2,
                                   const Double_t sn2)
{
  this->OscillateFD(dm2,sn2);
  //Get the true tau spectrum
  TH1D trueTauSpectrum = this->InverseOscillate(*fEffCorTau_FD,dm2,sn2);
  //true to reco the numus:
  this->TrueToRecoFD();
  //true to reco the taus:
  TH1D recoTauSpectrum = this->TrueToReco(trueTauSpectrum,
                                          *fRecoVsTrueEnergyTau_FD);
  //Add in the NCs
  this->PurityCorrectFDBackgroundNoOsc();
  //Add in the taus
  TH1D finalPrediction = this->AddInImpurity(*fRecoEnergyPred_FD,
                                             recoTauSpectrum);
  
  //This line is if I don't want taus:
  //       TH1D finalPrediction = *fRecoEnergyPred_FD;
  
  //Put the answer somewhere I can get to it later
  if (fRecoEnergyPred_FD){delete fRecoEnergyPred_FD; fRecoEnergyPred_FD=0;}
  fRecoEnergyPred_FD = new TH1D(finalPrediction);
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDDecaySpectrumBackgroundNoOsc(const Double_t alpha,
                                        const Double_t sn2)
{
  this->DecayFD(alpha,sn2);
  this->TrueToRecoFD();
  this->PurityCorrectFDBackgroundNoOsc();
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDDecoherenceSpectrumBackgroundNoOsc(const Double_t mu2,
                                              const Double_t sn2)
{
  this->DecohereFD(mu2,sn2);
  this->TrueToRecoFD();
  this->PurityCorrectFDBackgroundNoOsc();
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod::
PredictFDSpectrumBackgroundOscUncombined(const Double_t dm2,
                                         const Double_t sn2,
                                         const Double_t backdm2,
                                         const Double_t backsn2)
{
  this->OscillateFD(dm2,sn2);
  this->TrueToRecoFD();
  this->PurityCorrectFDBackgroundOscUncombined(backdm2,backsn2);
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDSpectrumAppearance(const Double_t dm2,
                              const Double_t sn2,
                              const Double_t backdm2,
                              const Double_t backsn2,
                              const Double_t appearanceFraction)
{
  this->OscillateFD(dm2,sn2);
  this->AddInAppearedEvents(backdm2,backsn2,appearanceFraction);
  this->TrueToRecoFD();
  this->PurityCorrectFDBackgroundOscCombined(backdm2,backsn2);
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
void NuMatrixMethod
::AddInAppearedEvents(const Double_t otherDm2,
                      const Double_t otherSn2,
                      const Double_t appearanceFraction)
{
  for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++) {
   
    Double_t energy = fTrueEnergyOsc_FD->GetBinCenter(i);
    Double_t transProb = otherSn2*(1-NuUtilities::OscillationWeight(energy, otherDm2, otherSn2));
    transProb *= appearanceFraction;
    
    fTrueEnergyOsc_FD->
      AddBinContent(i,
                    fAppearedTrueEffCor_FD->GetBinContent(i)*transProb);

    Double_t binError =
      TMath::Power(fTrueEnergyOsc_FD->GetBinError(i),2);
    binError += TMath::Power
      (fAppearedTrueEffCor_FD->GetBinContent(i)*transProb,2);
    if (binError>0){binError = TMath::Sqrt(binError);}
    else{binError = 0;}
    fTrueEnergyOsc_FD->SetBinError(i,binError);

  }
  return;
}


//____________________________________________________________________72
void NuMatrixMethod::SetFDAppearedFidFlux(const TH1D& hFDAppearedFidFlux)
{
  if (fAppearedTrueEffCor_FD){
    delete fAppearedTrueEffCor_FD;
    fAppearedTrueEffCor_FD = 0;
  }
  fAppearedTrueEffCor_FD = new TH1D(hFDAppearedFidFlux);
  fAppearedTrueEffCor_FD->SetDirectory(0);

  double val, err, energy;
  double newval, newerr;
  double eff, efferr, xsec;

  for(int i = 1; i <= fAppearedTrueEffCor_FD->GetNbinsX(); i++) {
    val = fAppearedTrueEffCor_FD->GetBinContent(i);
    err = fAppearedTrueEffCor_FD->GetBinError(i);
    energy = fAppearedTrueEffCor_FD->GetBinCenter(i);

    xsec = fXSec_CC_Graph->Eval(energy,0,"");
    eff = fEfficiency_FD->GetBinContent(i);
    efferr = fEfficiency_FD->GetBinError(i);
    
    // Add in xsec, eff Flux => Events
    newval = val * xsec * eff;
    
    newerr = TMath::Power(err/val,2) + TMath::Power(efferr/eff,2);
    newerr = TMath::Sqrt(newerr) * newval;
    
    fAppearedTrueEffCor_FD->SetBinContent(i,newval);
    fAppearedTrueEffCor_FD->SetBinError(i,newerr);
  }
  return;
}


//____________________________________________________________________72
const TH1D NuMatrixMethod::GetFDPotentialAppearanceFidFlux() const
{
  TH1D hPotential("hTrueEnergyPotentialAppearance",
                 "True energy potential appearance spectrum",
                 fNTrueBins,fTrueBinEdges);
  double val, err, energy;
  double newval, newerr;
  double eff, efferr, xsec;

  for(int i = 1; i <= hPotential.GetNbinsX(); i++) {
    val = fTrueEnergyEffCor_FD->GetBinContent(i);
    err = fTrueEnergyEffCor_FD->GetBinError(i);
    energy = hPotential.GetBinCenter(i);
    
    xsec = fXSec_CC_Graph->Eval(energy,0,"");
    eff = fEfficiency_FD->GetBinContent(i);
    efferr = fEfficiency_FD->GetBinError(i);
    
    // Remove xsec, eff: Events => Flux
    if (eff == 0 || xsec == 0) newval = 0;
    else newval = val / xsec / eff;
    
    newerr = TMath::Power(err/val,2) + TMath::Power(efferr/eff,2);
    newerr = TMath::Sqrt(newerr) * newval;
    
    hPotential.SetBinContent(i,newval);
    hPotential.SetBinError(i,newerr);
  }
  return hPotential;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDSpectrumBackgroundOscCombined(const Double_t dm2,
                                         const Double_t sn2,
                                         const Double_t backdm2,
                                         const Double_t backsn2)
{
  //applies osc to "T"
  this->OscillateFD(dm2,sn2);

  //does true to reco on osc true spectrum
  this->TrueToRecoFD();

  this->PurityCorrectFDBackgroundOscCombined(backdm2,backsn2);

  //give back the final oscillated reco energy spectrum with all
  //contaminations
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDDecaySpectrumBackgroundDecayCombined(const Double_t alpha,
                                                const Double_t sn2,
                                                const Double_t backalpha,
                                                const Double_t backsn2)
{
  //sn2 is now sin^{2}theta, not sin^{2}2theta

  //applies osc to "T"
  this->DecayFD(alpha,sn2);

  //does true to reco on osc true spectrum
  this->TrueToRecoFD();

  this->PurityCorrectFDBackgroundDecayCombined(backalpha,backsn2);

  //give back the final oscillated reco energy spectrum with all
  //contaminations
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D* NuMatrixMethod
::PredictFDDecoherenceSpectrumBackgroundDecohereCombined(const Double_t mu2,
                                                         const Double_t sn2,
                                                         const Double_t backmu2,
                                                         const Double_t backsn2)
{
  //sn2 is now sin^{2}theta, not sin^{2}2theta

  //applies osc to "T"
  this->DecohereFD(mu2,sn2);

  //does true to reco on osc true spectrum
  this->TrueToRecoFD();

  this->PurityCorrectFDBackgroundDecohereCombined(backmu2,backsn2);

  //give back the final oscillated reco energy spectrum with all
  //contaminations
  return fRecoEnergyPred_FD;
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::InverseOscillate(const TH1D& inputSpectrum,
                                            const Double_t dm2,
                                            const Double_t sn2) const
{
  TH1D outputSpectrum(inputSpectrum);
  outputSpectrum.Reset();

  for(int i=1;i<=outputSpectrum.GetNbinsX();i++) {
    Double_t energy = outputSpectrum.GetBinCenter(i);
    Double_t oscProb = 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
    outputSpectrum.
      SetBinContent(i,inputSpectrum.GetBinContent(i)*oscProb);
    outputSpectrum.
      SetBinError(i,inputSpectrum.GetBinError(i)*oscProb);
  }

  return outputSpectrum;
}

//____________________________________________________________________72
void NuMatrixMethod::OscillateFD(const Double_t dm2,
                                 const Double_t sn2)
{
  //8: Oscillate FD spectrum:
  for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++) {
    Double_t energy = fTrueEnergyOsc_FD->GetBinCenter(i);
    Double_t oscProb = NuUtilities::OscillationWeight(energy, dm2, sn2);
    fTrueEnergyOsc_FD->
      SetBinContent(i,fTrueEnergyEffCor_FD->GetBinContent(i)*oscProb);
    fTrueEnergyOsc_FD->
      SetBinError(i,fTrueEnergyEffCor_FD->GetBinError(i)*oscProb);
  }
}

// //____________________________________________________________________72
// void NuMatrixMethod::OscillateFD(const Double_t dm2,
//                               const Double_t sn2)
// {
//   //8: Oscillate FD spectrum:
//   for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++) {
//     Double_t energy = fTrueEnergyOsc_FD->GetBinCenter(i);
//     Double_t oscProb = NuUtilities::OscillationWeight(energy, dm2, sn2);
//     fTrueEnergyOsc_FD->
//       SetBinContent(i,fTrueEnergyEffCor_FD->GetBinContent(i)*oscProb);
//     fTrueEnergyOsc_FD->
//       SetBinError(i,fTrueEnergyEffCor_FD->GetBinError(i)*oscProb);
//   }
// }

//____________________________________________________________________72
void NuMatrixMethod::DecayFD(const Double_t alpha,
                             const Double_t sn2)
{
  //8: decay FD spectrum:
  for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++) {
    Double_t energy = fTrueEnergyOsc_FD->GetBinCenter(i);
    Double_t decayProb =
      TMath::Power
      (sn2 + (1.0-sn2)*exp(-(alpha*fFDDistance/Munits::km)/(2*energy)),
       2);
    fTrueEnergyOsc_FD->
      SetBinContent(i,fTrueEnergyEffCor_FD->GetBinContent(i)*decayProb);
    fTrueEnergyOsc_FD->
      SetBinError(i,fTrueEnergyEffCor_FD->GetBinError(i)*decayProb);
  }
}

//____________________________________________________________________72
void NuMatrixMethod::DecohereFD(const Double_t mu2,
                                const Double_t sn2)
{
  //8: Oscillate FD spectrum:
  for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++) {
    Double_t energy = fTrueEnergyOsc_FD->GetBinCenter(i);
    Double_t decoherenceProb =
      1.0 - (sn2/2)*
      (1 - exp(-(mu2*fFDDistance/Munits::km)/(2*energy)));
    fTrueEnergyOsc_FD->
      SetBinContent(i,fTrueEnergyEffCor_FD->GetBinContent(i)*decoherenceProb);
    fTrueEnergyOsc_FD->
      SetBinError(i,fTrueEnergyEffCor_FD->GetBinError(i)*decoherenceProb);
  }
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::TrueToReco(const TH1D& inputSpectrum,
                                      const TH2D& recoVsTrueHistogram) const
{
  TH1D outputSpectrum(inputSpectrum);
  outputSpectrum.Reset();

  //working space:
  Double_t *val = new Double_t[outputSpectrum.GetNbinsX()+1];
  Double_t *valerr = new Double_t[outputSpectrum.GetNbinsX()+1];
  
  //9: Convert to reco:
  for(int i=1;i<=outputSpectrum.GetNbinsX();i++){
    val[i-1] = 0; valerr[i-1] = 0;
  }
  for(int i=1;i<=inputSpectrum.GetNbinsX();i++){ //loop over true
    for(int j=1;j<=outputSpectrum.GetNbinsX()+1;j++){ //loop over reco
      val[j-1] += ( inputSpectrum.GetBinContent(i) *
                    recoVsTrueHistogram.GetBinContent(i,j));
      
      Double_t error = 0;
      if(inputSpectrum.GetBinContent(i)>0 &&
         recoVsTrueHistogram.GetBinContent(i,j)>0) {
        error = TMath::Power(inputSpectrum.GetBinError(i)/
                             inputSpectrum.GetBinContent(i),2);
        error += TMath::Power(recoVsTrueHistogram.GetBinError(i,j)/
                              recoVsTrueHistogram.GetBinContent(i,j),2);
        error *= TMath::Power(inputSpectrum.GetBinContent(i) * 
                              recoVsTrueHistogram.GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=outputSpectrum.GetNbinsX()+1;i++) {
    outputSpectrum.SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) outputSpectrum.SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else outputSpectrum.SetBinError(i,0);
  }

  return outputSpectrum;
}

//____________________________________________________________________72
void NuMatrixMethod::TrueToRecoFD()
{
  //working space:
  Double_t *val = new Double_t[fRecoEnergyOsc_FD->GetNbinsX()+1];
  Double_t *valerr = new Double_t[fRecoEnergyOsc_FD->GetNbinsX()+1];
  
  //9: Convert to reco:
  for(int i=1;i<=fRecoEnergyOsc_FD->GetNbinsX();i++){
    val[i-1] = 0; valerr[i-1] = 0;
  }
  for(int i=1;i<=fTrueEnergyOsc_FD->GetNbinsX();i++){ //loop over true
    for(int j=1;j<=fRecoEnergyOsc_FD->GetNbinsX()+1;j++){ //loop over reco
      val[j-1] += ( fTrueEnergyOsc_FD->GetBinContent(i) *
                    fRecoVsTrueEnergy_FD->GetBinContent(i,j));
      
      Double_t error = 0;
      if(fTrueEnergyOsc_FD->GetBinContent(i)>0 &&
         fRecoVsTrueEnergy_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(fTrueEnergyOsc_FD->GetBinError(i)/
                             fTrueEnergyOsc_FD->GetBinContent(i),2);
        error += TMath::Power(fRecoVsTrueEnergy_FD->GetBinError(i,j)/
                              fRecoVsTrueEnergy_FD->GetBinContent(i,j),2);
        error *= TMath::Power(fTrueEnergyOsc_FD->GetBinContent(i) * 
                              fRecoVsTrueEnergy_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=fRecoEnergyOsc_FD->GetNbinsX()+1;i++) {
    fRecoEnergyOsc_FD->SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) fRecoEnergyOsc_FD->SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else fRecoEnergyOsc_FD->SetBinError(i,0);
  }
}

//____________________________________________________________________72
void NuMatrixMethod::PurityCorrectFDBackgroundNoOsc()
{
  //10: Correct for purity (need to use noosc histogram to add 
  //    back in background otherwise it will be underestimated!):
  // THIS ONLY WORKS UNDER THE ASSUMPTION THAT YOUR BACKGROUND
  // DOES NOT OSCILLATE!
  // Background = ND, nutau CC, numubar CC, nue CC
  
  if (NuMMExtrapolation::kModularNoChargeCutFit == fextrpScheme ||
      NuMMExtrapolation::kModularPRLCCFit == fextrpScheme){
    for(int i=1;i<=fRecoEnergyPred_FD->GetNbinsX()+1;i++){
      fRecoEnergyPred_FD->
        SetBinContent(i,
                      fRecoEnergyOsc_FD->GetBinContent(i) +
                      fRecoUnoscNCBackground_FD->GetBinContent(i)
                      );
      Double_t error = 0.0;
      error +=
        TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
      error +=
        TMath::Power(fRecoUnoscNCBackground_FD->GetBinError(i),2);
      if (error > 0){
        fRecoEnergyPred_FD->SetBinError(i, TMath::Sqrt(error));
      }
      else{
        fRecoEnergyPred_FD->SetBinError(i, 0);
      }
    }
  }
  else{
    for(int i=1;i<=fRecoEnergyPred_FD->GetNbinsX()+1;i++){
      fRecoEnergyPred_FD->
        SetBinContent(i,fRecoEnergyOsc_FD->GetBinContent(i) +
                      fRecoEnergyNoOscPred_FD->GetBinContent(i) *
                      (1. - fPurity_FD->GetBinContent(i)));
      
      
      Double_t error = 0;
      if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
         (1.-fPurity_FD->GetBinContent(i))>0) {
        
        Double_t error =
          TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i)/
                       fRecoEnergyNoOscPred_FD->GetBinContent(i),
                       2);
        error += TMath::Power(fPurity_FD->GetBinError(i) /
                              (1. - fPurity_FD->GetBinContent(i)),
                              2);
        error = error*TMath::Power(fRecoEnergyNoOscPred_FD->GetBinContent(i) *
                                   (1. - fPurity_FD->GetBinContent(i)),
                                   2);
      }
      error += TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
      if(error>0) fRecoEnergyPred_FD->SetBinError(i,TMath::Sqrt(error));
      else fRecoEnergyPred_FD->SetBinError(i,0);   
    }
  }
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::AddInImpurity(const TH1D& inputSpectrum,
                                         const TH1D& impuritySpectrum) const
{
  TH1D outputSpectrum(inputSpectrum);
  outputSpectrum.Reset();
  
  for(int i=1;i<=outputSpectrum.GetNbinsX()+1;i++){
    outputSpectrum.
      SetBinContent(i,
                    inputSpectrum.GetBinContent(i) +
                    impuritySpectrum.GetBinContent(i)
                    );
    Double_t error = 0.0;
    error +=
      TMath::Power(inputSpectrum.GetBinError(i),2);
    error +=
      TMath::Power(impuritySpectrum.GetBinError(i),2);
    if (error > 0){
      outputSpectrum.SetBinError(i, TMath::Sqrt(error));
    }
    else{
      outputSpectrum.SetBinError(i, 0);
    }
  }
  
  return outputSpectrum;
}

//____________________________________________________________________72
void NuMatrixMethod::
PurityCorrectFDBackgroundOscUncombined(Double_t dm2Back,
                                        Double_t sn2Back)
{
  //10: Correct for purity assuming the background can oscillate.
  //Perform reco->true->oscillate->reco on the CC background.

  //working space:
  Double_t *val = new Double_t[fNTrueBins+1];
  Double_t *valerr = new Double_t[fNTrueBins+1];

  //A: Calculate unoscillated FD CC contamination:
  TH1D hRecoUnoscCCBackground_FD
    ("hUnoscCCBackground_FD",
     "Reco energy unoscillated CC background prediction",
     fNRecoBins,fRecoBinEdges);

  for(int i=1;i<=fNRecoBins+1;i++){
    hRecoUnoscCCBackground_FD.
      SetBinContent(i,
                    fRecoEnergyNoOscPred_FD->GetBinContent(i)*
                    fCCContamination_FD->GetBinContent(i)
                    );

    if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
       fCCContamination_FD->GetBinContent(i)>0) {
      Double_t error =
        TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i) / 
                     fRecoEnergyNoOscPred_FD->GetBinContent(i),2);
      error += TMath::Power(fCCContamination_FD->GetBinError(i) / 
                            fCCContamination_FD->GetBinContent(i),2);
      error =
        TMath::Sqrt(error) * hRecoUnoscCCBackground_FD.GetBinContent(i);
      hRecoUnoscCCBackground_FD.SetBinError(i,error);   
    }
    else hRecoUnoscCCBackground_FD.SetBinError(i,0);
  }
  
  //B: Get true energy spectrum of FD unoscillated CC contamination:
  TH1D hTrueUnoscCCBackground_FD
    ("fTrueUnoscCCBackground_FD",
     "True energy unoscillated CC background prediction",
     fNTrueBins,fTrueBinEdges);

  for(int i=1;i<=fNTrueBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNRecoBins+1;i++){
    //loop over reco dimension
    //include overflow
    for(int j=1;j<=fNTrueBins;j++){
      //loop over true dimension
      val[j-1] += ( hRecoUnoscCCBackground_FD.GetBinContent(i) * 
                    fRecoToTrueCCContamination_FD->GetBinContent(j,i) );
      Double_t error = 0;
      if(hRecoUnoscCCBackground_FD.GetBinContent(i)>0 &&
         fRecoToTrueCCContamination_FD->GetBinContent(j,i)>0) {
        error =
          TMath::Power(hRecoUnoscCCBackground_FD.GetBinError(i)/
                       hRecoUnoscCCBackground_FD.GetBinContent(i),2);
        error +=
          TMath::Power(fRecoToTrueCCContamination_FD->GetBinError(j,i)/
                       fRecoToTrueCCContamination_FD->GetBinContent(j,i),2);
        error *=
          TMath::Power(hRecoUnoscCCBackground_FD.GetBinContent(i) * 
                       fRecoToTrueCCContamination_FD->GetBinContent(j,i),2);
      }
      valerr[j-1] += error;
    }
  }
  for(int i=1;i<=fNTrueBins;i++) {
    hTrueUnoscCCBackground_FD.SetBinContent(i,val[i-1]);
    hTrueUnoscCCBackground_FD.SetBinError(i,TMath::Sqrt(valerr[i-1]));
  }

  //C: Oscillated CC background:
  TH1D hTrueOscCCBackground_FD
    ("hTrueOscCCBackground_FD",
     "True energy oscillated CC background prediction",
     fNTrueBins,fTrueBinEdges);
  for(int i=1;i<=fNTrueBins;i++) {
    Double_t energy = hTrueOscCCBackground_FD.GetBinCenter(i);
    Double_t oscProb = NuUtilities::OscillationWeight(energy, dm2Back, sn2Back);
    hTrueOscCCBackground_FD.
      SetBinContent(i,hTrueUnoscCCBackground_FD.GetBinContent(i)*oscProb);
    hTrueOscCCBackground_FD.
      SetBinError(i,hTrueUnoscCCBackground_FD.GetBinError(i)*oscProb);
  }

  //D: True->Reco CC background:
  TH1D hRecoOscCCBackground_FD
    ("hRecoOscCCBackground_FD",
     "Reco energy oscillated CC background prediction",
     fNRecoBins,fRecoBinEdges);
  
  for(int i=1;i<=fNRecoBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){ //loop over true
    for(int j=1;j<=fNRecoBins+1;j++){ //loop over reco
      val[j-1] += ( hTrueOscCCBackground_FD.GetBinContent(i) *
                    fTrueToRecoCCContamination_FD->GetBinContent(i,j));
      
      Double_t error = 0;
      if(hTrueOscCCBackground_FD.GetBinContent(i)>0 &&
         fTrueToRecoCCContamination_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(hTrueOscCCBackground_FD.GetBinError(i)/
                             hTrueOscCCBackground_FD.GetBinContent(i),2);
        error += TMath::Power(fTrueToRecoCCContamination_FD->GetBinError(i,j)/
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
        error *= TMath::Power(hTrueOscCCBackground_FD.GetBinContent(i) * 
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=fNRecoBins+1;i++) {
    hRecoOscCCBackground_FD.SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) hRecoOscCCBackground_FD.SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else hRecoOscCCBackground_FD.SetBinError(i,0);
  }

  //E: Calculate unoscillated FD NC contamination:
  /*
  TH1D hRecoUnoscNCBackground_FD
    ("hUnoscNCBackground_FD",
     "Reco energy unoscillated NC background prediction",
     fNRecoBins,fRecoBinEdges);

  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoUnoscNCBackground_FD->
      SetBinContent(i,
                    fRecoEnergyNoOscPred_FD->GetBinContent(i)*
                    fNCContamination_FD->GetBinContent(i)
                    );

    if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
       fNCContamination_FD->GetBinContent(i)>0) {
      Double_t error =
        TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i) / 
                     fRecoEnergyNoOscPred_FD->GetBinContent(i),2);
      error += TMath::Power(fNCContamination_FD->GetBinError(i) / 
                            fNCContamination_FD->GetBinContent(i),2);
      error =
        TMath::Sqrt(error) * fRecoUnoscNCBackground_FD->GetBinContent(i);
      fRecoUnoscNCBackground_FD->SetBinError(i,error);  
    }
    else fRecoUnoscNCBackground_FD->SetBinError(i,0);
  }
  */
  //F: Add the NC & oscillated CC backgrounds to the FD prediction:
  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoEnergyPred_FD->
      SetBinContent(i,
                    fRecoEnergyOsc_FD->GetBinContent(i) +
                    fRecoUnoscNCBackground_FD->GetBinContent(i) +
                    hRecoOscCCBackground_FD.GetBinContent(i)
                    );
    Double_t error = 0.0;
    error +=
      TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
    error +=
      TMath::Power(fRecoUnoscNCBackground_FD->GetBinError(i),2);
    error +=
      TMath::Power(hRecoOscCCBackground_FD.GetBinError(i),2);
    if (error > 0){
      fRecoEnergyPred_FD->SetBinError(i, TMath::Sqrt(error));
    }
    else{
      fRecoEnergyPred_FD->SetBinError(i, 0);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::
SetFDCCTrueUnoscBackground(const TH1D& hFDCCTrueUnoscBack)
{
  if (fSuppliedTrueUnoscCCBackground_FD){
    delete fSuppliedTrueUnoscCCBackground_FD;
    fSuppliedTrueUnoscCCBackground_FD = 0;
  }
  fSuppliedTrueUnoscCCBackground_FD = new TH1D(hFDCCTrueUnoscBack);
  fSuppliedTrueUnoscCCBackground_FD->SetDirectory(0);
  return;
}

//____________________________________________________________________72
void NuMatrixMethod::
PurityCorrectFDBackgroundOscCombined(Double_t dm2Back,
                                     Double_t sn2Back)
{
  //10: Correct for purity assuming the background can oscillate.
  //Use true CC background calculated from the extrapolation of the
  //other spectrum.

  //working space:
  Double_t *val = new Double_t[fNRecoBins+1];
  Double_t *valerr = new Double_t[fNRecoBins+1];

  //C: Oscillated CC background:
  TH1D hTrueOscCCBackground_FD
    ("hTrueOscCCBackground_FD",
     "True energy oscillated CC background prediction",
     fNTrueBins,fTrueBinEdges);
  for(int i=1;i<=fNTrueBins;i++) {
    Double_t energy = hTrueOscCCBackground_FD.GetBinCenter(i);
    Double_t oscProb = NuUtilities::OscillationWeight(energy, dm2Back, sn2Back);
    hTrueOscCCBackground_FD.
      SetBinContent(i,fSuppliedTrueUnoscCCBackground_FD->GetBinContent(i)*oscProb);
    hTrueOscCCBackground_FD.
      SetBinError(i,fSuppliedTrueUnoscCCBackground_FD->GetBinError(i)*oscProb);
  }

  //D: True->Reco CC background:
  TH1D hRecoOscCCBackground_FD
    ("hRecoOscCCBackground_FD",
     "Reco energy oscillated CC background prediction",
     fNRecoBins,fRecoBinEdges);
  
  for(int i=1;i<=fNRecoBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){ //loop over true
    for(int j=1;j<=fNRecoBins+1;j++){ //loop over reco
      val[j-1] += ( hTrueOscCCBackground_FD.GetBinContent(i) *
                    fTrueToRecoCCContamination_FD->GetBinContent(i,j));
      
      Double_t error = 0;
      if(hTrueOscCCBackground_FD.GetBinContent(i)>0 &&
         fTrueToRecoCCContamination_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(hTrueOscCCBackground_FD.GetBinError(i)/
                             hTrueOscCCBackground_FD.GetBinContent(i),2);
        error += TMath::Power(fTrueToRecoCCContamination_FD->GetBinError(i,j)/
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
        error *= TMath::Power(hTrueOscCCBackground_FD.GetBinContent(i) * 
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=fNRecoBins+1;i++) {
    hRecoOscCCBackground_FD.SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) hRecoOscCCBackground_FD.SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else hRecoOscCCBackground_FD.SetBinError(i,0);
  }

  //E: Calculate unoscillated FD NC contamination:
  /*
    //moved to CalculateNCContamination()
  TH1D hRecoUnoscNCBackground_FD
    ("hUnoscNCBackground_FD",
     "Reco energy unoscillated NC background prediction",
     fNRecoBins,fRecoBinEdges);

  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoUnoscNCBackground_FD->
      SetBinContent(i,
                    fRecoEnergyNoOscPred_FD->GetBinContent(i)*
                    fNCContamination_FD->GetBinContent(i)
                    );

    if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
       fNCContamination_FD->GetBinContent(i)>0) {
      Double_t error =
        TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i) / 
                     fRecoEnergyNoOscPred_FD->GetBinContent(i),2);
      error += TMath::Power(fNCContamination_FD->GetBinError(i) / 
                            fNCContamination_FD->GetBinContent(i),2);
      error =
        TMath::Sqrt(error) * fRecoUnoscNCBackground_FD->GetBinContent(i);
      fRecoUnoscNCBackground_FD->SetBinError(i,error);  
    }
    else fRecoUnoscNCBackground_FD->SetBinError(i,0);
  }
  */

  //F: Add the NC & oscillated CC backgrounds to the FD prediction:
  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoEnergyPred_FD->
      SetBinContent(i,
                    fRecoEnergyOsc_FD->GetBinContent(i) +
                    fRecoUnoscNCBackground_FD->GetBinContent(i) +
                    hRecoOscCCBackground_FD.GetBinContent(i)
                    );
    Double_t error = 0.0;
    error +=
      TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
    error +=
      TMath::Power(fRecoUnoscNCBackground_FD->GetBinError(i),2);
    error +=
      TMath::Power(hRecoOscCCBackground_FD.GetBinError(i),2);
    if (error > 0){
      fRecoEnergyPred_FD->SetBinError(i, TMath::Sqrt(error));
    }
    else{
      fRecoEnergyPred_FD->SetBinError(i, 0);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::
PurityCorrectFDBackgroundDecayCombined(Double_t alphaBack,
                                       Double_t sn2Back)
{
  //10: Correct for purity assuming the background can oscillate.
  //Use true CC background calculated from the extrapolation of the
  //other spectrum.

  //working space:
  Double_t *val = new Double_t[fNRecoBins+1];
  Double_t *valerr = new Double_t[fNRecoBins+1];

  //C: Oscillated CC background:
  TH1D hTrueOscCCBackground_FD
    ("hTrueOscCCBackground_FD",
     "True energy oscillated CC background prediction",
     fNTrueBins,fTrueBinEdges);
  for(int i=1;i<=fNTrueBins;i++) {
    Double_t energy = hTrueOscCCBackground_FD.GetBinCenter(i);
    Double_t oscProb =
      TMath::Power
      (sn2Back + (1.0-sn2Back)*exp(-(alphaBack*fFDDistance/Munits::km)/(2*energy)),
       2);
    hTrueOscCCBackground_FD.
      SetBinContent(i,fSuppliedTrueUnoscCCBackground_FD->GetBinContent(i)*oscProb);
    hTrueOscCCBackground_FD.
      SetBinError(i,fSuppliedTrueUnoscCCBackground_FD->GetBinError(i)*oscProb);
  }

  //D: True->Reco CC background:
  TH1D hRecoOscCCBackground_FD
    ("hRecoOscCCBackground_FD",
     "Reco energy oscillated CC background prediction",
     fNRecoBins,fRecoBinEdges);
  
  for(int i=1;i<=fNRecoBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){ //loop over true
    for(int j=1;j<=fNRecoBins+1;j++){ //loop over reco
      val[j-1] += ( hTrueOscCCBackground_FD.GetBinContent(i) *
                    fTrueToRecoCCContamination_FD->GetBinContent(i,j));
      
      Double_t error = 0;
      if(hTrueOscCCBackground_FD.GetBinContent(i)>0 &&
         fTrueToRecoCCContamination_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(hTrueOscCCBackground_FD.GetBinError(i)/
                             hTrueOscCCBackground_FD.GetBinContent(i),2);
        error += TMath::Power(fTrueToRecoCCContamination_FD->GetBinError(i,j)/
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
        error *= TMath::Power(hTrueOscCCBackground_FD.GetBinContent(i) * 
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=fNRecoBins+1;i++) {
    hRecoOscCCBackground_FD.SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) hRecoOscCCBackground_FD.SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else hRecoOscCCBackground_FD.SetBinError(i,0);
  }

  //F: Add the NC & oscillated CC backgrounds to the FD prediction:
  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoEnergyPred_FD->
      SetBinContent(i,
                    fRecoEnergyOsc_FD->GetBinContent(i) +
                    fRecoUnoscNCBackground_FD->GetBinContent(i) +
                    hRecoOscCCBackground_FD.GetBinContent(i)
                    );
    Double_t error = 0.0;
    error +=
      TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
    error +=
      TMath::Power(fRecoUnoscNCBackground_FD->GetBinError(i),2);
    error +=
      TMath::Power(hRecoOscCCBackground_FD.GetBinError(i),2);
    if (error > 0){
      fRecoEnergyPred_FD->SetBinError(i, TMath::Sqrt(error));
    }
    else{
      fRecoEnergyPred_FD->SetBinError(i, 0);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::
PurityCorrectFDBackgroundDecohereCombined(Double_t mu2Back,
                                          Double_t sn2Back)
{
  //10: Correct for purity assuming the background can oscillate.
  //Use true CC background calculated from the extrapolation of the
  //other spectrum.

  //working space:
  Double_t *val = new Double_t[fNRecoBins+1];
  Double_t *valerr = new Double_t[fNRecoBins+1];

  //C: Oscillated CC background:
  TH1D hTrueOscCCBackground_FD
    ("hTrueOscCCBackground_FD",
     "True energy oscillated CC background prediction",
     fNTrueBins,fTrueBinEdges);
  for(int i=1;i<=fNTrueBins;i++) {
    Double_t energy = hTrueOscCCBackground_FD.GetBinCenter(i);
    Double_t oscProb =
      1.0 - (sn2Back/2)*
      (1 - exp(-(mu2Back*fFDDistance/Munits::km)*(2*energy)));
    hTrueOscCCBackground_FD.
      SetBinContent(i,fSuppliedTrueUnoscCCBackground_FD->GetBinContent(i)*oscProb);
    hTrueOscCCBackground_FD.
      SetBinError(i,fSuppliedTrueUnoscCCBackground_FD->GetBinError(i)*oscProb);
  }

  //D: True->Reco CC background:
  TH1D hRecoOscCCBackground_FD
    ("hRecoOscCCBackground_FD",
     "Reco energy oscillated CC background prediction",
     fNRecoBins,fRecoBinEdges);
  
  for(int i=1;i<=fNRecoBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){ //loop over true
    for(int j=1;j<=fNRecoBins+1;j++){ //loop over reco
      val[j-1] += ( hTrueOscCCBackground_FD.GetBinContent(i) *
                    fTrueToRecoCCContamination_FD->GetBinContent(i,j));
      
      Double_t error = 0;
      if(hTrueOscCCBackground_FD.GetBinContent(i)>0 &&
         fTrueToRecoCCContamination_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(hTrueOscCCBackground_FD.GetBinError(i)/
                             hTrueOscCCBackground_FD.GetBinContent(i),2);
        error += TMath::Power(fTrueToRecoCCContamination_FD->GetBinError(i,j)/
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
        error *= TMath::Power(hTrueOscCCBackground_FD.GetBinContent(i) * 
                              fTrueToRecoCCContamination_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
      
    }
  }
  for(int i=1;i<=fNRecoBins+1;i++) {
    hRecoOscCCBackground_FD.SetBinContent(i,val[i-1]);
    if(valerr[i-1]>0) hRecoOscCCBackground_FD.SetBinError(i,TMath::Sqrt(valerr[i-1]));
    else hRecoOscCCBackground_FD.SetBinError(i,0);
  }

  //F: Add the NC & oscillated CC backgrounds to the FD prediction:
  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoEnergyPred_FD->
      SetBinContent(i,
                    fRecoEnergyOsc_FD->GetBinContent(i) +
                    fRecoUnoscNCBackground_FD->GetBinContent(i) +
                    hRecoOscCCBackground_FD.GetBinContent(i)
                    );
    Double_t error = 0.0;
    error +=
      TMath::Power(fRecoEnergyOsc_FD->GetBinError(i),2);
    error +=
      TMath::Power(fRecoUnoscNCBackground_FD->GetBinError(i),2);
    error +=
      TMath::Power(hRecoOscCCBackground_FD.GetBinError(i),2);
    if (error > 0){
      fRecoEnergyPred_FD->SetBinError(i, TMath::Sqrt(error));
    }
    else{
      fRecoEnergyPred_FD->SetBinError(i, 0);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::ResetUnoscHistograms()
{
  fRecoEnergyPurCor_ND->Reset();
  fTrueEnergyPurCor_ND->Reset();
  fTrueEnergyEffCor_ND->Reset();
  fTrueEnergyFlux_ND->Reset();
  fTrueEnergyMatrix_FD->Reset();
  fTrueEnergyCC_FD->Reset();
  fTrueEnergyEffCor_FD->Reset();
  fRecoEnergyNoOsc_FD->Reset();
  fRecoEnergyNoOscPred_FD->Reset();
  return;
}
 
//____________________________________________________________________72
void NuMatrixMethod::PurityCorrectND()
{ 
  //1: Correct ND reco spectrum for purity
  for(int i=1;i<=fNRecoBins+1;i++){
    fRecoEnergyPurCor_ND->
      SetBinContent(i,fRecoEnergyMeas_ND->GetBinContent(i) *
                    fPurity_ND->GetBinContent(i));
    if(fRecoEnergyMeas_ND->GetBinContent(i)>0 &&
       fPurity_ND->GetBinContent(i)>0) {
      Double_t error =
        TMath::Power(fRecoEnergyMeas_ND->GetBinError(i) / 
                     fRecoEnergyMeas_ND->GetBinContent(i),2);
      error += TMath::Power(fPurity_ND->GetBinError(i) / 
                            fPurity_ND->GetBinContent(i),2);
      error =
        TMath::Sqrt(error) * fRecoEnergyPurCor_ND->GetBinContent(i);
      fRecoEnergyPurCor_ND->SetBinError(i,error);       
    }
    else fRecoEnergyPurCor_ND->SetBinError(i,0);
  }
  return;
}

//____________________________________________________________________72
void NuMatrixMethod::RecoToTrueND()
{ 
  //working space:
  Double_t *val = new Double_t[fNTrueBins+1];
  Double_t *valerr = new Double_t[fNTrueBins+1];

  //2: Convert ND spectrum to true energies:
  for(int i=1;i<=fNTrueBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){
    //loop over reco dimension
    //include overflow
    for(int j=1;j<=fNRecoBins;j++){
      //loop over true dimension
      val[j-1] += ( fRecoEnergyPurCor_ND->GetBinContent(i) * 
                    fRecoVsTrueEnergy_ND->GetBinContent(j,i) );
      Double_t error = 0;
      if(fRecoEnergyPurCor_ND->GetBinContent(i)>0 &&
         fRecoVsTrueEnergy_ND->GetBinContent(j,i)>0) {
        error =
          TMath::Power(fRecoEnergyPurCor_ND->GetBinError(i)/
                       fRecoEnergyPurCor_ND->GetBinContent(i),2);
        error +=
          TMath::Power(fRecoVsTrueEnergy_ND->GetBinError(j,i)/
                       fRecoVsTrueEnergy_ND->GetBinContent(j,i),2);
        error *=
          TMath::Power(fRecoEnergyPurCor_ND->GetBinContent(i) * 
                       fRecoVsTrueEnergy_ND->GetBinContent(j,i),2);
      }
      valerr[j-1] += error;
    }
  }
  for(int i=1;i<=fNTrueBins;i++) {
    fTrueEnergyPurCor_ND->SetBinContent(i,val[i-1]);
    fTrueEnergyPurCor_ND->SetBinError(i,TMath::Sqrt(valerr[i-1]));
  }
  return;
}

//____________________________________________________________________72
void NuMatrixMethod::EfficiencyCorrectND()
{
  //3: Correct ND true energy spectrum for efficiencies:
  for(int i=1;i<=fNTrueBins;i++) {
    if(fEfficiency_ND->GetBinContent(i)>0 &&
       fTrueEnergyPurCor_ND->GetBinContent(i)>0) {
      fTrueEnergyEffCor_ND->
        SetBinContent(i,
                      fTrueEnergyPurCor_ND->GetBinContent(i) /
                      fEfficiency_ND->GetBinContent(i));        
      Double_t error =
        TMath::Power(fTrueEnergyPurCor_ND->GetBinError(i)/
                     fTrueEnergyPurCor_ND->GetBinContent(i),2);
      error += TMath::Power(fEfficiency_ND->GetBinError(i)/
                            fEfficiency_ND->GetBinContent(i),2);
      error =
        TMath::Sqrt(error) * fTrueEnergyEffCor_ND->GetBinContent(i);
      fTrueEnergyEffCor_ND->SetBinError(i,error);
    }
    else {
      fTrueEnergyEffCor_ND->SetBinContent(i,0);
      fTrueEnergyEffCor_ND->SetBinError(i,0);
    }
  }
  return;
}

//____________________________________________________________________72
void NuMatrixMethod::XSecCorrectND()
{
  //4: Get ND neutrino flux:
  if (fDoXSecStep){
    for(int i=1;i<=fNTrueBins;i++) {
      fTrueEnergyFlux_ND->
        SetBinContent(i,fTrueEnergyEffCor_ND->GetBinContent(i) *
                      (fFluxPoT/fNDPoT) / 
                      ( fNDFidMass *
                        fXSec_CC_Graph->Eval(fTrueEnergyFlux_ND->
                                             GetBinCenter(i),0,"") ) ); 
      fTrueEnergyFlux_ND->
        SetBinError(i,fTrueEnergyEffCor_ND->GetBinError(i) *
                    (fFluxPoT/fNDPoT) /
                    ( fNDFidMass *
                      fXSec_CC_Graph->Eval(fTrueEnergyFlux_ND->
                                           GetBinCenter(i),0,"") ) );
    }
  }
  else{
    for(int i=1;i<=fNTrueBins;i++) {
      fTrueEnergyFlux_ND->
        SetBinContent(i,fTrueEnergyEffCor_ND->GetBinContent(i) *
                      (fFluxPoT/fNDPoT) / 
                      ( fNDFidMass ) ); 
      fTrueEnergyFlux_ND->
        SetBinError(i,fTrueEnergyEffCor_ND->GetBinError(i) *
                    (fFluxPoT/fNDPoT) /
                    ( fNDFidMass ) );
    }
  }
  return;
}
 
//____________________________________________________________________72
void NuMatrixMethod::ExtrapolateNDToFD()
{
  //working space:
  Double_t *val = new Double_t[fNTrueBins+1];
  Double_t *valerr = new Double_t[fNTrueBins+1];

  //5: Get FD spectrum:
  for(int i=1;i<=fNTrueBins;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++) { //loop over ND bins
    for(int j=1;j<=fNTrueBins;j++) { //loop over FD bins
      val[j-1] += ( fTrueEnergyFlux_ND->GetBinContent(i) *
                    fFDVsNDMatrixRW->GetBinContent(i,j));
      Double_t error = 0;
      if(fTrueEnergyFlux_ND->GetBinContent(i)>0 &&
         fFDVsNDMatrixRW->GetBinContent(i,j)>0) {
        error = TMath::Power(fTrueEnergyFlux_ND->GetBinError(i)/
                             fTrueEnergyFlux_ND->GetBinContent(i),2);
        error += TMath::Power(fFDVsNDMatrixRW->GetBinError(i,j)/
                              fFDVsNDMatrixRW->GetBinContent(i,j),2);
        error *= TMath::Power(fTrueEnergyFlux_ND->GetBinContent(i) * 
                              fFDVsNDMatrixRW->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
    }
  }
  for(int i=1;i<=fNTrueBins;i++) {
    fTrueEnergyMatrix_FD->SetBinContent(i,val[i-1]);
    fTrueEnergyMatrix_FD->SetBinError(i,TMath::Sqrt(valerr[i-1]));
  }
  return;
}
//____________________________________________________________________72
void NuMatrixMethod::XSecCorrectFD()
{ 
  //6: Apply x-sec to get CC flux
  if (fDoXSecStep){
    for(int i=1;i<=fNTrueBins;i++) {
      fTrueEnergyCC_FD->SetBinContent(i,fTrueEnergyMatrix_FD->GetBinContent(i) *
                                      fXSec_CC_Graph->Eval(fTrueEnergyCC_FD->
                                                           GetBinCenter(i),0,"")*
                                      fFDFidMass * (fFDPoT/fFluxPoT));
      fTrueEnergyCC_FD->SetBinError(i,fTrueEnergyMatrix_FD->GetBinError(i) *
                                    fXSec_CC_Graph->Eval(fTrueEnergyCC_FD->
                                                         GetBinCenter(i),0,"")*
                                    fFDFidMass * (fFDPoT/fFluxPoT));
    }
  }
  else{
    for(int i=1;i<=fNTrueBins;i++) {
      fTrueEnergyCC_FD->SetBinContent(i,fTrueEnergyMatrix_FD->GetBinContent(i) *
                                      fFDFidMass * (fFDPoT/fFluxPoT));
      fTrueEnergyCC_FD->SetBinError(i,fTrueEnergyMatrix_FD->GetBinError(i) *
                                    fFDFidMass * (fFDPoT/fFluxPoT));
    }
  }
  return;
}

//____________________________________________________________________72
void NuMatrixMethod::PotentialTausTrueEnergy()
{
  TH1D hTrueEnergyTau_FD(this->TauXSecCorrectFD(*fTrueEnergyMatrix_FD));
  if (fEffCorTau_FD){delete fEffCorTau_FD; fEffCorTau_FD = 0;}
  fEffCorTau_FD = new TH1D(this->MultiplyInInefficiencies(hTrueEnergyTau_FD,
                                                          *fEfficiencyTau_FD));
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::TauXSecCorrectFD(const TH1D& tauFlux) const
{ 
  TH1D hTrueEnergyTau_FD("hTrueEnergyTau_FD","hTrueEnergyTau_FD",
                         fNTrueBins,fTrueBinEdges);
  //6: Apply x-sec to get CC flux
  if (fDoXSecStep){
    for(int i=1;i<=fNTrueBins;i++) {
      hTrueEnergyTau_FD.SetBinContent(i,tauFlux.GetBinContent(i) *
                                      fTau_XSec_CC_Graph->Eval(hTrueEnergyTau_FD.
                                                           GetBinCenter(i),0,"")*
                                      fFDFidMass * (fFDPoT/fFluxPoT));
      hTrueEnergyTau_FD.SetBinError(i,tauFlux.GetBinError(i) *
                                    fTau_XSec_CC_Graph->Eval(hTrueEnergyTau_FD.
                                                         GetBinCenter(i),0,"")*
                                    fFDFidMass * (fFDPoT/fFluxPoT));
    }
  }
  else{
    for(int i=1;i<=fNTrueBins;i++) {
      hTrueEnergyTau_FD.SetBinContent(i,tauFlux.GetBinContent(i) *
                                      fFDFidMass * (fFDPoT/fFluxPoT));
      hTrueEnergyTau_FD.SetBinError(i,tauFlux.GetBinError(i) *
                                    fFDFidMass * (fFDPoT/fFluxPoT));
    }
  }
  return hTrueEnergyTau_FD;
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::MultiplyInInefficiencies(const TH1D& inputSpectrum,
                                                    const TH1D& efficiencyPlot) const
{
  TH1D outputSpectrum(inputSpectrum);
  outputSpectrum.Reset();

  for(int i=1;i<=inputSpectrum.GetNbinsX();i++) {
    outputSpectrum.SetBinContent(i,inputSpectrum.GetBinContent(i) *
                                 efficiencyPlot.GetBinContent(i));
    if(inputSpectrum.GetBinContent(i)>0 &&
       efficiencyPlot.GetBinContent(i)>0) {
      Double_t error = TMath::Power(inputSpectrum.GetBinError(i)/
                                    inputSpectrum.GetBinContent(i),2);
      error += TMath::Power(efficiencyPlot.GetBinError(i)/
                            efficiencyPlot.GetBinContent(i),2);
      error = TMath::Sqrt(error)*outputSpectrum.GetBinContent(i);
      outputSpectrum.SetBinError(i,error);
    }
    else outputSpectrum.SetBinError(i,0);
  }
  return outputSpectrum;
}

//____________________________________________________________________72
void NuMatrixMethod::EfficiencyCorrectFD()
{ 
  //7: Apply FD efficiencies:
  for(int i=1;i<=fNTrueBins;i++) {
    fTrueEnergyEffCor_FD->SetBinContent(i,fTrueEnergyCC_FD->GetBinContent(i) *
                                        fEfficiency_FD->GetBinContent(i));
    if(fTrueEnergyCC_FD->GetBinContent(i)>0 &&
       fEfficiency_FD->GetBinContent(i)>0) {
      Double_t error = TMath::Power(fTrueEnergyCC_FD->GetBinError(i)/
                                    fTrueEnergyCC_FD->GetBinContent(i),2);
      error += TMath::Power(fEfficiency_FD->GetBinError(i)/
                            fEfficiency_FD->GetBinContent(i),2);
      error = TMath::Sqrt(error)*fTrueEnergyEffCor_FD->GetBinContent(i);
      fTrueEnergyEffCor_FD->SetBinError(i,error);
    }
    else fTrueEnergyEffCor_FD->SetBinError(i,0);
  }
  return;
}

//____________________________________________________________________72
const TH1D NuMatrixMethod::GetFDOtherCCTrueUnoscBackground()
{ 
  TH1D hTrueEnergyOtherCCTrueUnoscBack
    ("hTrueEnergyOtherCCTrueUnoscBack",
     "True energy unoscillated CC background for the other neutrinos",
     fNTrueBins,fTrueBinEdges);

  for(int i=1;i<=fNTrueBins;i++) {
    hTrueEnergyOtherCCTrueUnoscBack.
      SetBinContent(i,fTrueEnergyCC_FD->GetBinContent(i) *
                    fOtherEfficiency_FD->GetBinContent(i));
    if(fTrueEnergyCC_FD->GetBinContent(i)>0 &&
       fOtherEfficiency_FD->GetBinContent(i)>0) {
      Double_t error = TMath::Power(fTrueEnergyCC_FD->GetBinError(i)/
                                    fTrueEnergyCC_FD->GetBinContent(i),2);
      error += TMath::Power(fOtherEfficiency_FD->GetBinError(i)/
                            fOtherEfficiency_FD->GetBinContent(i),2);
      error = TMath::Sqrt(error)*hTrueEnergyOtherCCTrueUnoscBack.GetBinContent(i);
      hTrueEnergyOtherCCTrueUnoscBack.SetBinError(i,error);
    }
    else hTrueEnergyOtherCCTrueUnoscBack.SetBinError(i,0);
  }
  return hTrueEnergyOtherCCTrueUnoscBack;
}

//____________________________________________________________________72
void NuMatrixMethod::UnoscTrueToRecoFD()
{
  //working space:
  Double_t *val = new Double_t[fNRecoBins+1];
  Double_t *valerr = new Double_t[fNRecoBins+1];

  // Now we are going to digress in order to make the unoscillated FD spectrum:
  for(int i=1;i<=fNRecoBins+1;i++) { val[i-1] = 0; valerr[i-1] = 0; }
  for(int i=1;i<=fNTrueBins;i++){ //loop over true
    for(int j=1;j<=fNRecoBins+1;j++){ //loop over reco
      val[j-1] += ( fTrueEnergyEffCor_FD->GetBinContent(i) *
                    fRecoVsTrueEnergy_FD->GetBinContent(i,j) );
      Double_t error = 0;
      if(fTrueEnergyEffCor_FD->GetBinContent(i)>0 &&
         fRecoVsTrueEnergy_FD->GetBinContent(i,j)>0) {
        error = TMath::Power(fTrueEnergyEffCor_FD->GetBinError(i)/
                             fTrueEnergyEffCor_FD->GetBinContent(i),2);
        error += TMath::Power(fRecoVsTrueEnergy_FD->GetBinError(i,j)/
                              fRecoVsTrueEnergy_FD->GetBinContent(i,j),2);
        error *= TMath::Power(fTrueEnergyEffCor_FD->GetBinContent(i) * 
                              fRecoVsTrueEnergy_FD->GetBinContent(i,j),2);
      }
      valerr[j-1] += error;
    }
  }
  for(int i=1;i<=fNRecoBins+1;i++) {
    fRecoEnergyNoOsc_FD->SetBinContent(i,val[i-1]);
    fRecoEnergyNoOsc_FD->SetBinError(i,TMath::Sqrt(valerr[i-1]));
  }
  return;
}
 
//____________________________________________________________________72

void NuMatrixMethod::CalculateNCContamination()
{
  if (NuMMExtrapolation::kPRLCCAnalysis == fextrpScheme ||
      NuMMExtrapolation::kCCNuMuNoOscBackground == fextrpScheme){
    for(int i=1;i<=fNRecoBins+1;i++){
      fRecoUnoscNCBackground_FD->
      SetBinContent(i,
                    fRecoEnergyNoOscPred_FD->GetBinContent(i)*
                    (1.0-fPurity_FD->GetBinContent(i))
                    );
      
      if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
         1.0-fPurity_FD->GetBinContent(i)>0) {
        Double_t error =
          TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i) / 
                       fRecoEnergyNoOscPred_FD->GetBinContent(i),2);
        error += TMath::Power((fPurity_FD->GetBinError(i)) / 
                              (1.0-fPurity_FD->GetBinContent(i)),2);
        error =
          TMath::Sqrt(error)*fRecoUnoscNCBackground_FD->GetBinContent(i);
        fRecoUnoscNCBackground_FD->SetBinError(i,error);        
      }
      else fRecoUnoscNCBackground_FD->SetBinError(i,0);
    }
  }
  else{
    for(int i=1;i<=fNRecoBins+1;i++){
      fRecoUnoscNCBackground_FD->
      SetBinContent(i,
                    fRecoEnergyNoOscPred_FD->GetBinContent(i)*
                    fNCContamination_FD->GetBinContent(i)
                    );
      
      if(fRecoEnergyNoOscPred_FD->GetBinContent(i)>0 &&
         fNCContamination_FD->GetBinContent(i)>0) {
        Double_t error =
          TMath::Power(fRecoEnergyNoOscPred_FD->GetBinError(i) / 
                       fRecoEnergyNoOscPred_FD->GetBinContent(i),2);
        error += TMath::Power(fNCContamination_FD->GetBinError(i) / 
                              fNCContamination_FD->GetBinContent(i),2);
        error =
          TMath::Sqrt(error)*fRecoUnoscNCBackground_FD->GetBinContent(i);
        fRecoUnoscNCBackground_FD->SetBinError(i,error);        
      }
      else fRecoUnoscNCBackground_FD->SetBinError(i,0);
    }
  }
}

//____________________________________________________________________72
void NuMatrixMethod::UnoscPurityCorrectFD()
{ 
  // Correct for purity:
  for(int i=1;i<=fNRecoBins+1;i++){
    if(fRecoEnergyNoOsc_FD->GetBinContent(i)>1e-5 &&
       fPurity_FD->GetBinContent(i)>1e-5) {
      fRecoEnergyNoOscPred_FD->SetBinContent(i,fRecoEnergyNoOsc_FD->GetBinContent(i)/
                                             fPurity_FD->GetBinContent(i));
      Double_t error = TMath::Power(fRecoEnergyNoOsc_FD->GetBinError(i) / 
                                    fRecoEnergyNoOsc_FD->GetBinContent(i),2);
      error += TMath::Power(fPurity_FD->GetBinError(i) / 
                            fPurity_FD->GetBinContent(i),2);
      error = TMath::Sqrt(error) * fRecoEnergyNoOscPred_FD->GetBinContent(i);
      fRecoEnergyNoOscPred_FD->SetBinError(i,error);
    }
    else {
      fRecoEnergyNoOscPred_FD->SetBinContent(i,0);
      fRecoEnergyNoOscPred_FD->SetBinError(i,0);
    }
  }
  return;
}

//____________________________________________________________________72

---