NuMatrixFitter.cxx

Go to the documentation of this file.

// $Id: NuMatrixFitter.cxx,v 1.52 2009/09/26 00:44:48 nickd Exp $
//
// Performs Matrix Method analyses.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
// $Log: NuMatrixFitter.cxx,v $
// Revision 1.52  2009/09/26 00:44:48  nickd
// Removed casting pointer-to-int on cout statement
//
// This caused warning/errors on 64-bit. Besides, when pushing a pointer to
// ostream it writes it in hex, which is better!
//
// Revision 1.51  2008/04/07 12:59:51  evans
// Updating the PRL fit so it rebins to 0.5 GeV bins for the fit, and
// fits for normalisation.
//
// Revision 1.50  2008/03/28 22:10:53  ahimmel
//
//
// Got rid of compilation warning.
//
// Revision 1.49  2008/03/18 21:16:43  ahimmel
//
//
// Changes all in the Transition fitting.  Now chi2 curves are produced and filled in the Output object.  "experiments" is an argument like it is in NuTransitionFitterMinuit.
//
// Revision 1.48  2008/03/07 12:46:44  evans
// It seems I never committed the code for performing the multi-run fit
// for alternative models.
//
// Revision 1.47  2008/02/28 15:37:06  ahimmel
//
//
// Changed cout's to MSG()'s so they could be readily distinguished from debugging statements.
//
// Revision 1.46  2008/02/26 18:52:33  ahimmel
//
//
// Fixing memory leak in DoCPTFit by removing the chi2 histogram that does not get used.
//
// Revision 1.45  2008/02/20 16:07:33  evans
// The histogram binning is now configured from the .xml files.
//
// Revision 1.44  2008/02/19 19:16:38  evans
// Fixing memory management problems in the multi-run fits by holding the
// MMInput objects in a vector of pointers rather than a vector of
// MMInput objects.
//
// Revision 1.43  2008/02/18 13:48:07  rodriges
// Mini bug fix: if using scaled fd (fake) data, pass the scaled pots (not the
// overridden pots) to the matrix method to produce the unoscillated prediction.
//
// Also remove some now-irrelevant debug output
//
// Revision 1.42  2008/02/16 18:02:46  evans
// A multi-run PRL-style extrapolation and fit.
//
// Revision 1.41  2008/02/16 17:28:36  evans
// First attempt at getting the MultiRunCPT fitter in. It will
// compile. I've not checked it runs yet.
//
// Revision 1.40  2008/02/13 20:16:06  ahimmel
//
//
// Added TransitionFit code and got it up and working.
//
// Revision 1.39  2008/02/13 11:43:38  rodriges
// Yet more parameters passed via the XML:
//
// scaledFDPoTs is the target PoTs to scale the fake data to (for making sensitivity
// plots). To be used in a new script at the normalization stage in TheMethod.sh,
// which I'll commit shortly.
//
// overriddenFDPoTs, if set to a value >0, is a value to use for the far detector
// (fake) data PoTs instead of the value from hTotalPots. This is necessary for MDCs,
// where the PoT count from the file is wrong (it uses the default pots-per-file for
// MC, whereas the MDCs have a different exposure, eg 2.5e20).
//
// Revision 1.38  2008/02/12 14:43:25  evans
// When an alternative model is fitted, ensure that the output histograms
// in NuMatrixOutput correspond to that model.
//
// Revision 1.37  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.36  2008/02/09 14:30:59  rodriges
// Change the NoChargeCutFit method to use the new grid search parameters
//
// Revision 1.35  2008/02/08 19:20:42  rodriges
// Changes to allow setting the min, max and granularity for grid search. The
// make_xmls.pl currently sets the same things for the nu and bar variants, but they can
// be set differently in the XML file itself if necessary
//
// Revision 1.34  2008/02/08 15:49:09  evans
// The charge separated CC analysis can now be fitted with an energy cut
// on the high energy positive events.
//
// Revision 1.33  2008/01/28 10:20:54  hartnell
//
// Now passing the simFlag and using the enums directly:
//
// Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
//                                          SimFlag::kMC,
//                                          NuCuts::kJJE1);
//
// Revision 1.32  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.31  2008/01/20 01:19:14  evans
// In the new regime of systematically shifting the data instead of the
// MC, the NuXMLConfig is carried by the data file, not the helper
// file. So the fitter now looks for the NuXMLConfig in the ND data
// file. If it can't find it there, it looks in the helper file.
//
// Revision 1.30  2008/01/20 01:09:03  evans
// New functions to do the PRL-style CC analysis.
//
// Revision 1.29  2008/01/18 07:42:15  evans
// Lots of bug fixes to get the no charge cut CC fit working.
//
// Revision 1.28  2008/01/18 03:23:47  evans
// Updating the no charge cut fit to use the hTotalPot histograms and
// configurable fiducial masses.
//
// Revision 1.27  2008/01/16 19:12:06  evans
// Switching to using hTotalPot for PoT counting, and the new NuCuts
// fiducial mass function (changed in WriteHistosForFitter only).
//
// Revision 1.26  2008/01/09 19:55:30  ahimmel
//
//
// Updated the two matrix fitter methods to have shared methods that have the same input and output (CPTFit and DoCPTFit), or as well as they can match since the minuit method is doing a 4 parameter fit and the grid search is a 2 parameter fit.  They've also been set up so that the actual fitting procedure does not use NuMatrixMethod but instead used NuMatrixInput objects only do get the predicted FD spectrum.
//
// Revision 1.25  2008/01/09 15:34:49  evans
// A new function to do the standard CC analysis, including the
// antineutrinos, with a charge sign cut. It's the standard numubar
// 4-parameter CPT extrapolation and fit, with the neutrino and
// antineutrino parameters constrained to be identical.
//
// Revision 1.24  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.23  2007/12/24 15:55:21  evans
// Adding a function to do a transition analysis for multiple runs.
//
// Revision 1.22  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.21  2007/12/19 15:32:53  evans
// Code to do the extrapolation and fitting separately for an arbitrary
// number of datasets (i.e. runs with different SKZP weights). Use
// WriteMultiRunHistosForFitter(...) and DoMultiRunCPTFit(...). There's a
// corresponding new constructor. Instead of passing in filenames, simply
// pass in vectors of the filenames for each run (in the same order for
// each type of file).
//
// Revision 1.20  2007/12/14 23:57:10  ahimmel
//
//
// Lots of small changes to make these files work well together in the automated matrix method.  Mostly changing the types of strings being passed and changing object references to object pointers and the like.  As much as possible I've tried to keep the interfaces back compatible.  The only times this won't work is when I've had to change the format in which an object is troed (reference to pointer for the FDData in NuMatrixOutput for example).
//
// Revision 1.19  2007/12/06 10:00:38  hartnell
//
// 2 fixes to make it compile.
//
// Firstly, comment these lines out:
//
//   //xmlConfig->SetDirectory(0);
//   //xmlConfig.Close();
//
// Secondly, fix a typo:
//
//   xmlCofig needs an "n".
//
// Revision 1.18  2007/12/06 01:13:53  ahimmel
//
//
// Pass the NuXMlConfig object from file to file.
//
// Revision 1.17  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.16  2007/12/02 07:49:53  rhatcher
// Add '#include <cmath>' so that round() is available on gcc 4.0.1.
//
// Revision 1.15  2007/11/30 19:41:06  ahimmel
//
//
// Adapted new XMLConfig-sensitive fitting routine into the framework (using DoCPTFit) so that code is not repeated.
//
// Revision 1.14  2007/11/30 16:40:08  evans
// Adding code to nu a numu->numubar appearance analysis with the matrix
// method.
//
// Revision 1.13  2007/11/29 22:25:56  ahimmel
//
//
// Made some modifications to work with the systematic study framework.
//
// NuMatrixOutput:
// added data members for the best fit dm2bar and sn2bar
//
// NuMatrixFitter:
 // Create and alternate inferface to CPTFit called DoCPTFit that uses NuMatrixInput and NuMatrixOutpout objects directly rather than filenames.
//
// NuXMLConfig:
// Added a FullName() function that returns a name that includes the shift value.
//
// Revision 1.12  2007/11/28 10:56:02  evans
// Adding a new version of the modular fitting function that takes a
// NuXMLConfig object instead of a string labelling the systematic shift.
//
// Revision 1.11  2007/11/20 16:36:03  evans
// Updating the NuMatrixOutput object to hold a no oscillation FD prediction.
//
// Revision 1.10  2007/11/20 15:57:09  evans
// Fixing the grid search as it was missing out the bins at the dm2 maximum.
//
// Revision 1.9  2007/11/20 12:41:18  evans
// Debugging the modular MM fitter.
//
// Revision 1.8  2007/11/19 16:27:43  evans
// Writing a function to do a nubar CPT modular MM fit for systematics
// studies.
//
// Revision 1.7  2007/11/16 18:31:52  hartnell
//
// Allow output file to be specified when writing histograms for the fitter
//
// Revision 1.6  2007/11/15 15:37:48  evans
// Making the grid searches a little more configurable. Making a chi2 surface.
//
// Revision 1.5  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.4  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.3  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.2  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.1  2007/09/06 12:23:23  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
//

#include <cmath>

#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TROOT.h"

#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuCuts.h"
#include "NtupleUtils/NuFluctuator.h"
#include "NtupleUtils/NuMatrixFitter.h"
#include "NtupleUtils/NuMatrixInput.h"
#include "NtupleUtils/NuMatrixMethod.h"
#include "NtupleUtils/NuMatrixOutput.h"
#include "NtupleUtils/NuUtilities.h"
#include "NtupleUtils/NuXMLConfig.h"
#include "Conventions/SimFlag.h"

ClassImp(NuMatrixFitter)

CVSID("$Id: NuMatrixFitter.cxx,v 1.52 2009/09/26 00:44:48 nickd Exp $");

//____________________________________________________________________72
NuMatrixFitter::NuMatrixFitter()
{
  fndDataFilename = "";
  ffdDataFilename = "";
  fhelperFilename = "";
  fxsecFilename = "";

  SetGridParamsDefault();
  fxmlConfig=0;
}

//____________________________________________________________________72
NuMatrixFitter::NuMatrixFitter(const Double_t FluxPoT,
                               const string ndData,
                               const string fdData,
                               const string helperFile,
                               const string xsecFile,
                               const NuXMLConfig* xmlConfig)
{
  //This constructor is the setup for extrapolating single runs
  ffluxPoT = FluxPoT;
  fndDataFilename = ndData;
  ffdDataFilename = fdData;
  fhelperFilename = helperFile;
  fxsecFilename = xsecFile;
  foutputFilename =
    "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutput.root";
  fxmlConfig=xmlConfig;
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
  } else {
    SetGridParamsDefault();
  }
}

//____________________________________________________________________72
NuMatrixFitter::NuMatrixFitter(const Double_t fluxPoT,
                               const vector<string> ndDataFiles,
                               const vector<string> fdDataFiles,
                               const vector<string> helperFiles,
                               const string xsecFile,
                               const NuXMLConfig* xmlConfig)
{
  //This constructor is the setup for extrapolating multiple runs. It
  //will work just as well on one run: just put one filename in each
  //vector.

  ffluxPoT = fluxPoT;
  fvndDataFilenames = ndDataFiles;
  fvfdDataFilenames = fdDataFiles;
  fvhelperFilenames = helperFiles;
  fxsecFilename = xsecFile;
  foutputFilename =
    "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutput.root";
  fxmlConfig=xmlConfig;
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
  } else {
    SetGridParamsDefault();
  }
}

//____________________________________________________________________72
NuMatrixFitter::~NuMatrixFitter()
{
}

//____________________________________________________________________72
void NuMatrixFitter::SetXMLConfig(const NuXMLConfig* xmlConfig) 
{
  fxmlConfig=xmlConfig;
  this->FillGridSearchParams(xmlConfig);
}

//____________________________________________________________________72
void NuMatrixFitter::SetGridParamsDefault()
{
  fdm2NuLow=0.001;
  fdm2NuHigh=0.01;
  fdm2NuGranularity=0.001;

  fsn2NuLow=0;
  fsn2NuHigh=1;
  fsn2NuGranularity=0.05;

  fdm2BarLow=0.001;
  fdm2BarHigh=0.01;
  fdm2BarGranularity=0.001;

  fsn2BarLow=0;
  fsn2BarHigh=1;
  fsn2BarGranularity=0.05;

  ftransitionProbLow=0;
  ftransitionProbHigh=1;
  ftransitionProbGranularity=0.05;

}
//____________________________________________________________________72
void NuMatrixFitter::FillGridSearchParams(const NuXMLConfig* xmlConfig)
{
   fdm2NuLow = xmlConfig->DM2NuLow();
   fdm2NuHigh = xmlConfig->DM2NuHigh();
   fdm2NuGranularity = xmlConfig->DM2NuGranularity();

   fsn2NuLow = xmlConfig->SN2NuLow();
   fsn2NuHigh = xmlConfig->SN2NuHigh();
   //MSG("NuMatrixFitter",Msg::kInfo)  << "SN2NuGranularity = " << xmlConfig->SN2NuGranularity() << endl;
   fsn2NuGranularity = xmlConfig->SN2NuGranularity();

   fdm2BarLow = xmlConfig->DM2BarLow();
   fdm2BarHigh = xmlConfig->DM2BarHigh();
   fdm2BarGranularity = xmlConfig->DM2BarGranularity();

   fsn2BarLow = xmlConfig->SN2BarLow();
   fsn2BarHigh = xmlConfig->SN2BarHigh();
   fsn2BarGranularity = xmlConfig->SN2BarGranularity();

   ftransitionProbLow = xmlConfig->TransitionProbLow();
   ftransitionProbHigh = xmlConfig->TransitionProbHigh();
   ftransitionProbGranularity = xmlConfig->TransitionProbGranularity();
}

//____________________________________________________________________72
const Double_t NuMatrixFitter::CalculateChi2(const TH1D* fdPred, 
                                                                                                                                                                                 const TH1D* fdData) const
{
  //Return Chi2.

  Double_t chi2 = 0;

  for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fdPred->GetBinContent(i);
    Double_t dnu = fdData->GetBinContent(i);
//     MSG("NuMatrixFitter",Msg::kInfo)  << "Bin " << i
//       << ": MC " << mnu
//       << "; data " << dnu
//       << endl;
    if (mnu<0.0001){mnu = 0.0001;}
    chi2 += (dnu-mnu)*(dnu-mnu)/mnu;
  }
  return chi2;
}

//____________________________________________________________________72
const Double_t NuMatrixFitter::CalculateLikelihood(const TH1D* fdPred,
                                                   const TH1D* fdData) const
{
  //Aim to minimise -2lnL. That's what I'm returning.

  Double_t like = 0;

  for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fdPred->GetBinContent(i);
    Double_t dnu = fdData->GetBinContent(i);
    if (mnu<0.0001){mnu = 0.0001;}
    if (dnu){
      like += 2*(mnu - dnu + dnu*log(dnu/mnu));
    }
    else{like += 2*(mnu - dnu);}
  }
//   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
  return like;
}

//____________________________________________________________________72
const Double_t NuMatrixFitter::CalculateLikelihoodNormPenalty(const TH1D* fdPred,
                                                              const TH1D* fdData,
                                                              const Double_t normalisation) const
{
  //Aim to minimise -2lnL. That's what I'm returning.

  Double_t normOneSigma = 0.04;
  Double_t like = 0;

  for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fdPred->GetBinContent(i);
    Double_t dnu = fdData->GetBinContent(i);
    if (mnu<0.0001){mnu = 0.0001;}
    if (dnu){
      like += 2*(mnu - dnu + dnu*log(dnu/mnu)) +
        ((normalisation - 1.0)*(normalisation-1.0))/
        (normOneSigma*normOneSigma);
    }
    else{like += 2*(mnu - dnu);}
  }
//   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
  return like;
}

//____________________________________________________________________72
const Double_t NuMatrixFitter::NormalisationPenaltyTerm(const Double_t normalisation) const
{
  Double_t normOneSigma = 0.04;
  return ((normalisation - 1.0)*(normalisation-1.0))/
    (normOneSigma*normOneSigma);
}

//____________________________________________________________________72
const Double_t NuMatrixFitter
::CalculateLikelihood(const TH1D* fdPred,
                      const TH1D* fdData,
                      const Double_t Ecut) const
{
  //Aim to minimise -2lnL. That's what I'm returning.

  Double_t like = 0;
  Int_t i = 1;
  if (i>0.0){
    i = fdPred->GetXaxis()->FindBin(Ecut);
  }
  for (; i<=fdPred->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fdPred->GetBinContent(i);
    Double_t dnu = fdData->GetBinContent(i);
    if (mnu<0.0001){mnu = 0.0001;}
    if (dnu){
      like += 2*(mnu - dnu + dnu*log(dnu/mnu));
    }
    else{like += 2*(mnu - dnu);}
  }
  return like;
}

//____________________________________________________________________72
void NuMatrixFitter::CCAnalysis()
{
  //Get the FD data & PoT histograms
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  //  TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyAll_FD");
  fdData->SetDirectory(0);
  fdData->Sumw2();
  TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
  hfdPoTTorTGT->SetDirectory(0);
  TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
  hfdSpillsPerFile->SetDirectory(0);
  TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
  hfdRun->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
  hndPoTTorTGT->SetDirectory(0);
  ndDataFile.Close();

  //Calculate the total PoT
  Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
   Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
   //   Double_t fdPoT = 2.5e20;
  MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;

  //Hard coded, unchecked fiducial masses. Should make configurable.
  Double_t ndFidMass = 45.128*Munits::tonne;
  //  Double_t fdFidMass = ndFidMass*(448/68)*(14.0-0.3*0.3/1.0);
  Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
  // Double_t fdFidMass = 5.404*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);

  NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        foutputFilename);
  mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNuMuNoOscBackground);
  //  mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNoChargeCut);
  mmNuMu.PredictFDFluxUnosc();
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  Double_t sn2 = 0.0;
  Double_t dm2 = 0.0;

  Int_t numSn2bins = (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity);
  Int_t numDm2bins = (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity);
  TH2D hChi2("hChi2",
             "",
             numSn2bins,
             fsn2NuLow,
             fsn2NuHigh,
             numDm2bins,
             fdm2NuLow,
             fdm2NuHigh);
   for (/*Double_t*/ sn2 = fsn2NuLow; sn2 < fsn2NuHigh+fsn2NuGranularity; sn2 += fsn2NuGranularity){
     for (/*Double_t*/ dm2 = fdm2NuLow; dm2 < fdm2NuHigh+fdm2NuGranularity; dm2 += fdm2NuGranularity){
      const TH1D* fdPred =
        mmNuMu.PredictFDSpectrumBackgroundNoOsc(dm2, sn2);
      Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2+fsn2NuGranularity/2.0);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2+fdm2NuGranularity/2.0);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      hChi2.SetBinContent(bin2D,chi2);
       MSG("NuMatrixFitter",Msg::kInfo)  << "Result: " << hChi2.GetBinContent(bin2D) << endl;
       MSG("NuMatrixFitter",Msg::kInfo)  << "dm2: " << dm2
           << "; sn2: " << sn2
           << "; chi2: " << chi2 << endl;
      if (chi2 < bestChi2){
        bestChi2 = chi2;
        bestsn2 = sn2;
        bestdm2 = dm2;
      }//if
     }//dm2
   }//sn2
  MSG("NuMatrixFitter",Msg::kInfo)  << "Best fit: dm2 = " << bestdm2
       << "; sn2 = " << bestsn2
       << "; bestChi2 = " << bestChi2 << endl;
  
  TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutput.root","RECREATE");
  fdData->Write();
  hChi2.Write();
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}
  
}

//____________________________________________________________________72
NuMatrixOutput* NuMatrixFitter::DoTransitionFit(//const NuMatrixInput mmInput,
                                                const char * inputFilename,
                                                const TH1D* numubarfdData,
                                                const Double_t dm2,
                                                const Double_t sn2) const
{
    static int calls = 0;
    TString hN = "";
    hN += calls;
    ++calls;
    NuMatrixInput *mmInput = new NuMatrixInput(inputFilename);
    
    Double_t bestChi2 = 1.0e10;
    Double_t bestTransitionProb = 0.0;
    
    float transitionProbLow = -0.5;
    float transitionProbHigh = 0.5;
    float transitionProbGranularity = 0.05;

    TH1D* numubarfdPred = 0;
    int bins = (int)((transitionProbHigh - transitionProbLow)/transitionProbGranularity+1);
    TH1D* surface = new TH1D("hLikelihood"+hN,"Likelihood Curve", bins, transitionProbLow, transitionProbHigh);
    int ibin = 1;
    for (Double_t transitionProb = transitionProbLow;
         transitionProb <= transitionProbHigh+transitionProbGranularity/2.0;
         transitionProb += transitionProbGranularity){
        numubarfdPred = mmInput->PredictFDSpectrumTransition(dm2, sn2, transitionProb);
        Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData);
        surface->SetBinContent(ibin++, chi2);

        
        if (chi2 < bestChi2){
            bestChi2 = chi2;
            bestTransitionProb = transitionProb;
        }
        delete numubarfdPred; numubarfdPred = 0;
    }//transitionProb
    
    MSG("NuMatrixFitter.cxx",Msg::kInfo) << "Best fit: transitionProb = " << bestTransitionProb << endl;
    
    NuMatrixOutput *mmOutput = new NuMatrixOutput();
    mmOutput->NuMuBarFDData(numubarfdData);
    
    mmOutput->NuMuBarBestFitFDPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, bestTransitionProb));
    mmOutput->NuMuBarNoTransPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, 0.0)); 
    mmOutput->NuMuBarNoOscPrediction(*mmInput->PredictFDSpectrumTransition(0.0, 0.0, 0.0));   
    mmOutput->NuMuChi2TransSurface(*surface);
    mmOutput->BestFitPoint(sn2, dm2, bestTransitionProb);
    
    delete mmInput;
    delete surface;
    return mmOutput;
}

//____________________________________________________________________72
NuMatrixOutput* NuMatrixFitter::DoCPTFit(NuMatrixInput *mmInput,
                                         TH1D* numubarfdData,
                                         const Double_t dm2,
                                         const Double_t sn2)
{
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2bar = 0.0;
  Double_t bestdm2bar = 0.0;
  
  //Create the Chi2 histogram
  Int_t numSn2BarBins = (Int_t) round((fsn2BarHigh-fsn2BarLow)/fsn2BarGranularity) + 1;
  Int_t numDm2BarBins = (Int_t) round((fdm2BarHigh-fdm2BarLow)/fdm2BarGranularity) + 1;
  TH2D hChi2Bar("hChi2Bar",
                "",
                numSn2BarBins,
                fsn2BarLow-fsn2BarGranularity/2.0,
                fsn2BarHigh+fsn2BarGranularity/2.0,
                numDm2BarBins,
                fdm2BarLow-fdm2BarGranularity/2.0,
                fdm2BarHigh+fdm2BarGranularity/2.0); 
  
  for (Double_t sn2bar = fsn2BarLow;
       sn2bar <= fsn2BarHigh+fsn2BarGranularity/2.0;
       sn2bar += fsn2BarGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "sn2bar = " << sn2bar << "; bestChi2 " << bestChi2 << endl;
    for (Double_t dm2bar = fdm2BarLow;
         dm2bar <= fdm2BarHigh+fdm2BarGranularity/2.0;
         dm2bar += fdm2BarGranularity){
      
      TH1D* numubarfdPred = mmInput->PredictFDSpectrumNuMuBar
        (dm2, sn2, dm2bar, sn2bar);
      
      Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData);
      
      //       chi2 += this->CalculateLikelihood(numufdPred,numufdData);
      Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
      Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
      Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
      hChi2Bar.SetBinContent(bin2D,chi2);
      if (chi2 < bestChi2){
        bestChi2 = chi2;
        bestsn2bar = sn2bar;
        bestdm2bar = dm2bar;
      }//if
    }//dm2bar
  }//sn2bar
  
  MSG("NuMatrixFitter.cxx",Msg::kInfo)
    << "Best fit:" << endl
    << "    sin2bar = " << bestsn2bar << endl
    << "    dm2bar = " << bestdm2bar << endl;
  
  NuMatrixOutput *mmOutput = new NuMatrixOutput();
  mmOutput->NuMuBarChi2Surface(hChi2Bar);
  mmOutput->NuMuBarFDData(numubarfdData);
  
  mmOutput->NuMuBarBestFitFDPrediction
    (*mmInput->PredictFDSpectrumNuMuBar(dm2, sn2, bestdm2bar, bestsn2bar));
  mmOutput->NuMuBarNoOscPrediction
    (*mmInput->PredictFDSpectrumNuMuBar(0., 0., 0., 0.));
  
  mmOutput->BestFitPoint(bestsn2bar, bestdm2bar);
  
        return mmOutput;
}

//____________________________________________________________________72
NuMatrixOutput* NuMatrixFitter
::DoCCFitChargeCut(const NuMatrixInput& mmInput,
                   const TH1D* numufdData,
                   const TH1D* numubarfdData,
                   const Double_t nubarEcut,
                   const NuModel_t model)
{
  NuMatrixMethod mmNuMu(mmInput,
                        NuMMExtrapolation::kModularNuMuCPTFit);
  NuMatrixMethod mmNuMuBar(mmInput,
                           NuMMExtrapolation::kModularNuMuBarCPTFit);
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  
  //Create the Chi2 histogram
  Int_t numSn2Bins =
    (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
  Int_t numDm2Bins =
    (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
  TH2D hChi2("hChi2",
             "",
             numSn2Bins,
             fsn2NuLow-fsn2NuGranularity/2.0,
             fsn2NuHigh+fsn2NuGranularity/2.0,
             numDm2Bins,
             fdm2NuLow-fdm2NuGranularity/2.0,
             fdm2NuHigh+fdm2NuGranularity/2.0); 
  
  for (Double_t sn2 = fsn2NuLow;
       sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
       sn2 += fsn2NuGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "sn2 = " << sn2 << "; best chi2 " << bestChi2 << endl;
    for (Double_t dm2 = fdm2NuLow;
         dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
         dm2 += fdm2NuGranularity){
      const TH1D* numufdPred = 0;
      const TH1D* numubarfdPred = 0;
      if (kOscillation == model){
        numufdPred =
          mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                        sn2,
                                                        dm2,
                                                        sn2);
        numubarfdPred =
          mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                           sn2,
                                                           dm2,
                                                           sn2);
      }
      else if (kDecay == model){
        numufdPred =
          mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(dm2,
                                                               sn2,
                                                               dm2,
                                                               sn2);
        numubarfdPred =
          mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(dm2,
                                                                  sn2,
                                                                  dm2,
                                                                  sn2);
//      MSG("NuMatrixFitter",Msg::kInfo)  << "Decayed; dm2 " << dm2
//           << ", sn2 " << sn2
//           << ", bin 4 contents: " << numufdPred->GetBinContent(4)
//           << endl;
      }
      else if (kDecoherence == model){
        numufdPred =
          mmNuMu.PredictFDDecoherenceSpectrumBackgroundDecohereCombined(dm2,
                                                                        sn2,
                                                                        dm2,
                                                                        sn2);
        numubarfdPred =
          mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundDecohereCombined(dm2,
                                                                           sn2,
                                                                           dm2,
                                                                           sn2);
      }
      else {
        MSG("NuMatrixFitter.cxx",Msg::kError)
          << "Unknown neutrino disappearance model!" << endl;
      }
      
      Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData,nubarEcut);
      chi2 += this->CalculateLikelihood(numufdPred,numufdData);
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      hChi2.SetBinContent(bin2D,chi2);
      if (chi2 < bestChi2){
        bestChi2 = chi2;
        bestsn2 = sn2;
        bestdm2 = dm2;
      }//if
    }//dm2
  }//sn2
  
  MSG("NuMatrixFitter.cxx",Msg::kInfo)
    << "Best fit:" << endl
    << "    sin2 = " << bestsn2 << endl
    << "    dm2 = " << bestdm2 << endl;
  
  NuMatrixOutput *mmOutput = new NuMatrixOutput();

  mmOutput->NuMuBarChi2Surface(hChi2);
  mmOutput->NuMuBarFDData(numubarfdData);
  if (kOscillation == model){  
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(bestdm2,
                                                         bestsn2,
                                                         bestdm2,
                                                         bestsn2));
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(0.0,
                                                         0.0,
                                                         0.0,
                                                         0.0));
    mmOutput->NuMuBestFitFDPrediction
      (*mmNuMu.PredictFDSpectrumBackgroundOscCombined(bestdm2,
                                                      bestsn2,
                                                      bestdm2,
                                                      bestsn2));
    mmOutput->NuMuNoOscPrediction
      (*mmNuMu.PredictFDSpectrumBackgroundOscCombined(0.0,
                                                      0.0,
                                                      0.0,
                                                      0.0));
  }
  else if (kDecay == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(bestdm2,
                                                                bestsn2,
                                                                bestdm2,
                                                                bestsn2));
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(0.0,
                                                                0.0,
                                                                0.0,
                                                                0.0));
    mmOutput->NuMuBestFitFDPrediction
      (*mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(bestdm2,
                                                             bestsn2,
                                                             bestdm2,
                                                             bestsn2));
    mmOutput->NuMuNoOscPrediction
      (*mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(0.0,
                                                             0.0,
                                                             0.0,
                                                             0.0));
  }
  else if (kDecoherence == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.
       PredictFDDecoherenceSpectrumBackgroundDecohereCombined(bestdm2,
                                                              bestsn2,
                                                              bestdm2,
                                                              bestsn2));
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.
       PredictFDDecoherenceSpectrumBackgroundDecohereCombined(0.0,
                                                              0.0,
                                                              0.0,
                                                              0.0));
    mmOutput->NuMuBestFitFDPrediction
      (*mmNuMu.
       PredictFDDecoherenceSpectrumBackgroundDecohereCombined(bestdm2,
                                                              bestsn2,
                                                              bestdm2,
                                                              bestsn2));
    mmOutput->NuMuNoOscPrediction
      (*mmNuMu.
       PredictFDDecoherenceSpectrumBackgroundDecohereCombined(0.0,
                                                              0.0,
                                                              0.0,
                                                              0.0));
  }
  else{
    MSG("NuMatrixFitter.cxx",Msg::kError)
      << "Unknown neutrino disappearance model!" << endl;
  }
  mmOutput->BestFitPoint(bestsn2, bestdm2);

  mmOutput->NuMuChi2Surface(hChi2);
  mmOutput->NuMuFDData(numufdData);  
  mmOutput->BestFitPoint(bestsn2, bestdm2);
  
  return mmOutput;
}

//____________________________________________________________________72
NuMatrixOutput* NuMatrixFitter::DoNoChargeCutFit(NuMatrixInput& mmInput,
                                                 TH1D* fdData,
                                                 const NuModel_t model)
{
  NuMatrixMethod mmNuMuBar(mmInput,
                           NuMMExtrapolation::kModularNoChargeCutFit);
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  
  //Create the Chi2 histogram
  Int_t numSn2Bins =
    (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
  Int_t numDm2Bins =
    (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
  TH2D hChi2("hChi2",
                "",
                numSn2Bins,
                fsn2NuLow-fsn2NuGranularity/2.0,
                fsn2NuHigh+fsn2NuGranularity/2.0,
                numDm2Bins,
                fdm2NuLow-fdm2NuGranularity/2.0,
                fdm2NuHigh+fdm2NuGranularity/2.0); 

  for (Double_t sn2 = fsn2NuLow;
       sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
       sn2 += fsn2NuGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "sn2 = " << sn2 << "; bestChi2 " << bestChi2 
      << " (bestdm2 " << bestdm2
      << ", bestsn2 " << bestsn2 << ")" << endl;
    for (Double_t dm2 = fdm2NuLow;
                                 dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
                                 dm2 += fdm2NuGranularity){
      const TH1D* fdPred = 0;
      if (kOscillation == model){
        fdPred =
          mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(dm2,
                                                     sn2);
      }
      else if (kDecay == model){
        fdPred =
          mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(dm2,
                                                          sn2);
      }
      else if (kDecoherence == model){
        fdPred =
          mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
                                                                sn2);
      }
      else {
        MSG("NuMatrixFitter.cxx",Msg::kInfo)
          << "Invalid neutrino disappearance model" << endl;
      }
      
      Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
      //       chi2 += this->CalculateLikelihood(numufdPred,numufdData);
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      hChi2.SetBinContent(bin2D,chi2);
      if (chi2 < bestChi2){
                                bestChi2 = chi2;
                                bestsn2 = sn2;
                                bestdm2 = dm2;
      }//if
    }//dm2
  }//sn2

  MSG("NuMatrixFitter.cxx",Msg::kInfo)
                << "Best fit:" << endl
                << "    sin2 = " << bestsn2 << endl
    << "    dm2 = " << bestdm2 << endl;
        
  NuMatrixOutput *mmOutput = new NuMatrixOutput();
  mmOutput->NuMuBarChi2Surface(hChi2);
  mmOutput->NuMuBarFDData(fdData);
  if (kOscillation == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(bestdm2,
                                                   bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(0.0,
                                                   0.0));
  }
  else if (kDecay == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
                                                        bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(0.0,
                                                        0.0));
  }
  else if (kDecoherence == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
                                                        bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
                                                              0.0));
  }
  else{
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "Invalid neutrino disappearance model" << endl;
  }
  mmOutput->BestFitPoint(bestsn2, bestdm2);
  
        return mmOutput;
}

//____________________________________________________________________72
NuMatrixOutput* NuMatrixFitter::DoPRLCCFit(const NuMatrixInput& mmInput,
                                           const TH1D* fdData,
                                           const NuModel_t model)
{
  NuMatrixMethod mmNuMuBar(mmInput,
                           NuMMExtrapolation::kModularPRLCCFit);
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  
  //Create the Chi2 histogram
  Int_t numSn2Bins =
    (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
  Int_t numDm2Bins =
    (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
  TH2D hChi2("hChi2",
                "",
                numSn2Bins,
                fsn2NuLow-fsn2NuGranularity/2.0,
                fsn2NuHigh+fsn2NuGranularity/2.0,
                numDm2Bins,
                fdm2NuLow-fdm2NuGranularity/2.0,
                fdm2NuHigh+fdm2NuGranularity/2.0); 

  for (Double_t sn2 = fsn2NuLow;
       sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
       sn2 += fsn2NuGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "sn2 = " << sn2 << "; bestChi2 " << bestChi2 
      << " (bestdm2 " << bestdm2
      << ", bestsn2 " << bestsn2 << ")" << endl;
    for (Double_t dm2 = fdm2NuLow;
         dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
         dm2 += fdm2NuGranularity){
      const TH1D* fdPred = 0;
      if (kOscillation == model){
        fdPred =
          mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(dm2,
                                                     sn2);
      }
      else if (kDecay == model){
        fdPred =
          mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(dm2,
                                                          sn2);
      }
      else if (kDecoherence == model){
        fdPred =
          mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
                                                                sn2);
      }
      else {
        MSG("NuMatrixFitter.cxx",Msg::kError)
          << "Invalid neutrino disappearance model" << endl;
      }
      // 
      //TH1D* outputhist = (TH1D*)fdPred->Clone();
      //outputhist->SetName(Form("fdpred_s%.2f_m%.2f", sn2, dm2*1000));
      //outputhist->Write();
      Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      hChi2.SetBinContent(bin2D,chi2);
      //MSG("NuMatrixFitter.cxx",Msg::kInfo)
      //  << "sn2=" << sn2 << " dm2=" << dm2 << " chi2=" << chi2 << endl;
      if (chi2 < bestChi2){
                                bestChi2 = chi2;
                                bestsn2 = sn2;
                                bestdm2 = dm2;
      }//if
    }//dm2
  }//sn2

  MSG("NuMatrixFitter.cxx",Msg::kInfo)
                << "Best fit:" << endl
                << "    sin2 = " << bestsn2 << endl
    << "    dm2 = " << bestdm2 << endl;
        
  NuMatrixOutput *mmOutput = new NuMatrixOutput();
  mmOutput->NuMuBarChi2Surface(hChi2);
  mmOutput->NuMuBarFDData(fdData);
  
  if (kOscillation == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(bestdm2,
                                                   bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(0.0,
                                                   0.0));
  }
  else if (kDecay == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
                                                        bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(0.0,
                                                        0.0));
  }
  else if (kDecoherence == model){
    mmOutput->NuMuBarBestFitFDPrediction
      (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
                                                        bestsn2));
    
    mmOutput->NuMuBarNoOscPrediction
      (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
                                                              0.0));
  }
  else{
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "Invalid neutrino disappearance model" << endl;
  }

  mmOutput->BestFitPoint(bestsn2, bestdm2);
  MSG("NuMatrixFitter.cxx",Msg::kInfo) << "returning mmoutput" << endl;
  return mmOutput;
}

//____________________________________________________________________72
const NuMatrixOutput* NuMatrixFitter
::DoMultiRunTransitionFit(const vector<NuMatrixInput> mmInput,
                          const vector<TH1D> numubarfdData,
                          const Double_t dm2,
                          const Double_t sn2) const
{
  vector<NuMatrixMethod*> nubarExtrapolators;
  for (vector<NuMatrixInput>::const_iterator inputIt = mmInput.begin();
       inputIt != mmInput.end();
       ++inputIt){
    nubarExtrapolators.push_back
      (new NuMatrixMethod(*inputIt,
                          NuMMExtrapolation::kModularNuMuBarTransitionFit));
  }

  Double_t bestChi2 = 1.0e10;
  Double_t bestTransitionProb = 0.0;

//   Double_t transitionProbLow = 0.0;
//   Double_t transitionProbHigh = 1.0;
//   Double_t transitionProbGranularity = 0.05;

  //Create the chi2 histogram
  Int_t numTransitionProbBins =
    (Int_t) round((ftransitionProbHigh-ftransitionProbLow)/
                  ftransitionProbGranularity) + 1;
  TH1D hChi2Trans("hChi2Trans",
                  "",
                  numTransitionProbBins,
                  ftransitionProbLow-ftransitionProbGranularity/2.0,
                  ftransitionProbHigh+ftransitionProbGranularity/2.0);
  for (Double_t transitionProb = ftransitionProbLow;
       transitionProb <= ftransitionProbHigh+ftransitionProbGranularity/2.0;
       transitionProb += ftransitionProbGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "Fitting transitionProb = " << transitionProb << endl;

      Double_t chi2 = 0.0;
      vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
      vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
      for (;
           mmNuBarIt != nubarExtrapolators.end();
           ++mmNuBarIt, ++fdDataIt){
        const TH1D* numubarfdPred =
          (*mmNuBarIt)->PredictFDSpectrumAppearance(dm2,
                                                    sn2,
                                                    dm2,
                                                    sn2,
                                                    transitionProb);  
        chi2 += this->CalculateLikelihood(numubarfdPred,&(*fdDataIt));
      }

    Int_t xBin = hChi2Trans.GetXaxis()->FindBin(transitionProb);
    hChi2Trans.SetBinContent(xBin,chi2);
    if (chi2 < bestChi2){
      bestChi2 = chi2;
      bestTransitionProb = transitionProb;
    }
  }//transitionProb

  MSG("NuMatrixFitter.cxx",Msg::kInfo)
    << "Best fit: transitionProb = " << bestTransitionProb << endl;
  
  TH1D* hbestFitPred = 0;
  TH1D* hnoTransPred = 0;
  TH1D* hnoOscPred = 0;
  TH1D* hcombinedFDData = 0;
  Bool_t firstTime = true;
  vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
  vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
  for (;
       mmNuBarIt != nubarExtrapolators.end();
       ++mmNuBarIt, ++fdDataIt){
    if (firstTime){
      hbestFitPred = new TH1D
        (*(*mmNuBarIt)->
         PredictFDSpectrumAppearance(dm2,
                                     sn2,
                                     dm2,
                                     sn2,
                                     bestTransitionProb));
      hnoTransPred = new TH1D
        (*(*mmNuBarIt)->
         PredictFDSpectrumAppearance(dm2,
                                     sn2,
                                     dm2,
                                     sn2,
                                     0.0));
      hnoOscPred = new TH1D
        (*(*mmNuBarIt)->
         PredictFDSpectrumAppearance(0.0,
                                     0.0,
                                     0.0,
                                     0.0,
                                     0.0));
      hcombinedFDData = new TH1D(*fdDataIt);
      firstTime = false;
    }
    else{
      hbestFitPred->Add
        ((*mmNuBarIt)->
         PredictFDSpectrumAppearance(dm2,
                                     sn2,
                                     dm2,
                                     sn2,
                                     bestTransitionProb));
      hnoTransPred->Add
        ((*mmNuBarIt)->
         PredictFDSpectrumAppearance(dm2,
                                     sn2,
                                     dm2,
                                     sn2,
                                     0.0));
      hnoOscPred->Add
        ((*mmNuBarIt)->
         PredictFDSpectrumAppearance(0.0,
                                     0.0,
                                     0.0,
                                     0.0,
                                     0.0));
      hcombinedFDData->Add(&(*fdDataIt));
    }
  }
        
  NuMatrixOutput *mmOutput = new NuMatrixOutput();
  mmOutput->NuMuBarChi2TransSurface(hChi2Trans);
  mmOutput->NuMuBarFDData(*hcombinedFDData);
  mmOutput->NuMuBarBestFitFDPrediction(*hbestFitPred);
  mmOutput->NuMuBarNoTransPrediction(*hnoTransPred);
  mmOutput->NuMuBarNoOscPrediction(*hnoOscPred);
  mmOutput->BestFitPoint(sn2, dm2, bestTransitionProb);
  return mmOutput;
}

//____________________________________________________________________72
void NuMatrixFitter::MultiRunPRLCCFit(const vector<string> inputFilenames,
                                      const vector<string> fdDataFilenames,
                                      const vector<string> outputFilenames,
                                      const string combinedOutputFilename,
                                      const NuModel_t model)
{
  MSG("NuMatrixFitter",Msg::kInfo)  << "NuMatrixFitter::PRLCCFit()" << endl;
  vector<TH1D> nuMuFDDataHistos;
  for (vector<string>::const_iterator fdDataIt = fdDataFilenames.begin();
       fdDataIt != fdDataFilenames.end();
       ++fdDataIt){
    TFile fdDataFile((*fdDataIt).c_str(),"READ");
    TH1D* nubarfdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
    nubarfdData->SetDirectory(0);
    nubarfdData->Sumw2();
    fdDataFile.Close();
    nuMuFDDataHistos.push_back(*nubarfdData);
  }
  
  vector<NuMatrixInput*> mmInputs;
  vector<NuXMLConfig> xmlConfigs;
  for (vector<string>::const_iterator inFileIt = inputFilenames.begin();
       inFileIt != inputFilenames.end();
       ++inFileIt){

    TFile infile((*inFileIt).c_str(), "READ");
    
    NuMatrixInput *mmInput = 0;
    mmInput = (NuMatrixInput*) infile.Get("NuMatrixInput");
    if (!mmInput) {
      mmInput = new NuMatrixInput(&infile);
    }
    mmInputs.push_back(mmInput);
    if (!fxmlConfig){
      fxmlConfig = (NuXMLConfig*) infile.Get("NuXMLConfig");
    }
    if (fxmlConfig) {
      FillGridSearchParams(fxmlConfig);
      MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
      xmlConfigs.push_back(*fxmlConfig);
    }
    else {
      MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
      NuXMLConfig dummy;
      xmlConfigs.push_back(dummy);
      //objName = "Nominal";
    }
  }
  
  pair<vector<NuMatrixOutput*>,NuMatrixOutput*> mmOutputs =
    this->DoMultiRunPRLCCFit(mmInputs,
                             nuMuFDDataHistos,
                             model);
   
  vector<string>::const_iterator outFileIt = outputFilenames.begin();
  vector<NuXMLConfig>::const_iterator xmlConfigIt = xmlConfigs.begin();
  for (vector<NuMatrixOutput*>::iterator mmOutIt = mmOutputs.first.begin();
       mmOutIt != mmOutputs.first.end();
       ++mmOutIt, ++outFileIt, ++xmlConfigIt){
    MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << *outFileIt << endl;
    TFile *outputFile = new TFile((*outFileIt).c_str(),"RECREATE");
    outputFile->cd();
    (*mmOutIt)->XMLConfig(&(*xmlConfigIt));
    (*mmOutIt)->Write("NuMatrixOutput");
    outputFile->Close();
  }
  TFile combinedOutFile(combinedOutputFilename.c_str(),"RECREATE");
  MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << combinedOutFile.GetName() << endl;
  combinedOutFile.cd();
  mmOutputs.second->XMLConfig(&(*xmlConfigs.begin()));
  mmOutputs.second->Write("NuMatrixOutput");
  combinedOutFile.Close();

     /*
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
  */
}

//____________________________________________________________________72
void NuMatrixFitter::MultiRunCPTFit(const vector<string> inputFilenames,
                                    const vector<string> fdDataFilenames,
                                    const Double_t dm2,
                                    const Double_t sn2,
                                    const vector<string> outputFilenames,
                                    const string combinedOutputFilename)
{
  vector<TH1D> nuMuBarFDDataHistos;
  for (vector<string>::const_iterator fdDataIt = fdDataFilenames.begin();
       fdDataIt != fdDataFilenames.end();
       ++fdDataIt){
    TFile fdDataFile((*fdDataIt).c_str(),"READ");
    TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
    numubarfdData->SetDirectory(0);
    numubarfdData->Sumw2();
    fdDataFile.Close();
    nuMuBarFDDataHistos.push_back(*numubarfdData);
  }
  
  vector<NuMatrixInput*> mmInputs;
  vector<NuXMLConfig> xmlConfigs;
  for (vector<string>::const_iterator inFileIt = inputFilenames.begin();
       inFileIt != inputFilenames.end();
       ++inFileIt){

    TFile infile((*inFileIt).c_str(), "READ");
    
    NuMatrixInput *mmInput = 0;
    mmInput = (NuMatrixInput*) infile.Get("NuMatrixInput");
    if (!mmInput) {
      mmInput = new NuMatrixInput(&infile);
    }
    mmInputs.push_back(mmInput);
    if (!fxmlConfig){
      fxmlConfig = (NuXMLConfig*) infile.Get("NuXMLConfig");
    }
    if (fxmlConfig) {
      FillGridSearchParams(fxmlConfig);
      MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
      xmlConfigs.push_back(*fxmlConfig);
    }
    else {
      MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
      NuXMLConfig dummy;
      xmlConfigs.push_back(dummy);
      //objName = "Nominal";
    }
  }
  
  pair<vector<NuMatrixOutput*>,NuMatrixOutput*> mmOutputs =
    this->DoMultiRunCPTFit(mmInputs,
                           nuMuBarFDDataHistos,
                           dm2,
                           sn2);
   
  vector<string>::const_iterator outFileIt = outputFilenames.begin();
  vector<NuXMLConfig>::const_iterator xmlConfigIt = xmlConfigs.begin();
  for (vector<NuMatrixOutput*>::iterator mmOutIt = mmOutputs.first.begin();
       mmOutIt != mmOutputs.first.end();
       ++mmOutIt, ++outFileIt, ++xmlConfigIt){
    MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << *outFileIt << endl;
    TFile *outputFile = new TFile((*outFileIt).c_str(),"RECREATE");
    outputFile->cd();
    (*mmOutIt)->XMLConfig(&(*xmlConfigIt));
    (*mmOutIt)->Write("NuMatrixOutput");
    outputFile->Close();
  }
  TFile combinedOutFile(combinedOutputFilename.c_str(),"RECREATE");
  combinedOutFile.cd();
  mmOutputs.second->XMLConfig(&(*xmlConfigs.begin()));
  mmOutputs.second->Write("NuMatrixOutput");
  combinedOutFile.Close();

     /*
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
  */
}

//____________________________________________________________________72
const pair<vector<NuMatrixOutput*>,NuMatrixOutput*>
NuMatrixFitter::DoMultiRunCPTFit(const vector<NuMatrixInput*> mmInput,
                                 const vector<TH1D> numubarfdData,
                                 const Double_t dm2,
                                 const Double_t sn2) const
{
  vector<NuMatrixMethod*> nubarExtrapolators;
  for (vector<NuMatrixInput*>::const_iterator inputIt = mmInput.begin();
       inputIt != mmInput.end();
       ++inputIt){
    nubarExtrapolators.push_back
      (new NuMatrixMethod(**inputIt,
                          NuMMExtrapolation::kModularNuMuBarCPTFit));
  }
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2bar = 0.0;
  Double_t bestdm2bar = 0.0;
  
  //Create the Chi2 histogram
  Int_t numSn2BarBins =
    (Int_t) round((fsn2BarHigh-fsn2BarLow)/fsn2BarGranularity) + 1;
  Int_t numDm2BarBins =
    (Int_t) round((fdm2BarHigh-fdm2BarLow)/fdm2BarGranularity) + 1;
  TH2D hChi2Bar("hChi2Bar",
                "",
                numSn2BarBins,
                fsn2BarLow-fsn2BarGranularity/2.0,
                fsn2BarHigh+fsn2BarGranularity/2.0,
                numDm2BarBins,
                fdm2BarLow-fdm2BarGranularity/2.0,
                fdm2BarHigh+fdm2BarGranularity/2.0); 
  
  for (Double_t sn2bar = fsn2BarLow;
       sn2bar <= fsn2BarHigh+fsn2BarGranularity/2.0;
       sn2bar += fsn2BarGranularity){
    MSG("NuMatrixFitter.cxx",Msg::kInfo)
      << "sn2bar = " << sn2bar << "; bestChi2 " << endl;
    for (Double_t dm2bar = fdm2BarLow;
         dm2bar <= fdm2BarHigh+fdm2BarGranularity/2.0;
         dm2bar += fdm2BarGranularity){

      Double_t chi2 = 0.0;
      vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
      vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
      for (;
           mmNuBarIt != nubarExtrapolators.end();
           ++mmNuBarIt, ++fdDataIt){
        const TH1D* numubarfdPred =
          (*mmNuBarIt)->PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                               sn2bar,
                                                               dm2,
                                                               sn2);
        chi2 += this->CalculateLikelihood(numubarfdPred,&(*fdDataIt));
      }
      
      Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
      Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
      Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
      hChi2Bar.SetBinContent(bin2D,chi2);
      if (chi2 < bestChi2){
        bestChi2 = chi2;
        bestsn2bar = sn2bar;
        bestdm2bar = dm2bar;
      }//if
    }//dm2bar
  }//sn2bar

  MSG("NuMatrixFitter.cxx",Msg::kInfo)
                << "Best fit:" << endl
                << "    sin2bar = " << bestsn2bar << endl
    << "    dm2bar = " << bestdm2bar << endl;
  
  TH1D* hbestFitPred = 0;
  TH1D* hnoOscPred = 0;
  TH1D* hcombinedFDData = 0;
  vector<NuMatrixOutput*> mmOuts;
  Bool_t firstTime = true;
  vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
  vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
  for (;
       mmNuBarIt != nubarExtrapolators.end();
       ++mmNuBarIt, ++fdDataIt){
    NuMatrixOutput* mmOutput = new NuMatrixOutput();
    mmOutput->NuMuBarChi2Surface(hChi2Bar);
    mmOutput->NuMuBarFDData(*fdDataIt);
    mmOutput->NuMuBarBestFitFDPrediction(*(*mmNuBarIt)->
                                         PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
                                                                                bestsn2bar,
                                                                                dm2,
                                                                                sn2));
    mmOutput->NuMuBarNoOscPrediction(*(*mmNuBarIt)->
                                     PredictFDSpectrumBackgroundOscCombined(0.0,
                                                                            0.0,
                                                                            0.0,
                                                                            0.0));
    mmOutput->BestFitPoint(bestsn2bar, bestdm2bar);
    mmOuts.push_back(mmOutput);
    if (firstTime){
      hbestFitPred = new TH1D
        (*(*mmNuBarIt)->
         PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
                                                bestsn2bar,
                                                dm2,
                                                sn2));
      hnoOscPred = new TH1D
        (*(*mmNuBarIt)->
         PredictFDSpectrumBackgroundOscCombined(0.0,
                                                0.0,
                                                0.0,
                                                0.0));
      hcombinedFDData = new TH1D(*fdDataIt);
      firstTime = false;
    }
    else{
      hbestFitPred->Add
        ((*mmNuBarIt)->
         PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
                                                bestsn2bar,
                                                dm2,
                                                sn2));
      hnoOscPred->Add
        ((*mmNuBarIt)->
         PredictFDSpectrumBackgroundOscCombined(0.0,
                                                0.0,
                                                0.0,
                                                0.0));
      hcombinedFDData->Add(&(*fdDataIt));
    }
  }
        
  NuMatrixOutput *mmOutputCombined = new NuMatrixOutput();
  mmOutputCombined->NuMuBarChi2Surface(hChi2Bar);
  mmOutputCombined->NuMuBarFDData(*hcombinedFDData);
  mmOutputCombined->NuMuBarBestFitFDPrediction(*hbestFitPred);
  mmOutputCombined->NuMuBarNoOscPrediction(*hnoOscPred);
  mmOutputCombined->BestFitPoint(bestsn2bar, bestdm2bar);

  pair<vector<NuMatrixOutput*>,NuMatrixOutput*> allMMOuts(mmOuts,mmOutputCombined);
  return allMMOuts;
}

//____________________________________________________________________72
TH1D
NuMatrixFitter::RebinForFit(const TH1D* fineHist)
{
//   NuBinningScheme::NuBinningScheme_t binningScheme = NuBinningScheme::k0_5GeVTo200;
//   const NuUtilities utils;
//   const vector<Double_t> recoBins = utils.RecoBins(binningScheme);
//   const Int_t numFitBins = recoBins.size()-1;
//   Double_t *fitBinsArray;
//   fitBinsArray=new Double_t[numFitBins+1];
//   {
//     Int_t i=0;
//     for (vector<Double_t>::const_iterator itBin = recoBins.begin();
//       itBin != recoBins.end();
//       ++itBin, ++i){
//       fitBinsArray[i] = *itBin;
//     }
//   }

  static Int_t stupidRoot = 0;
  TH1D newHist(Form("NewHist%d",stupidRoot++),"",800,0,200);
  newHist.SetBinContent(1,0);
  for (Int_t oldBin = 1; oldBin<fineHist->GetNbinsX()+1; ++oldBin){
    newHist.SetBinContent(oldBin+1,fineHist->GetBinContent(oldBin));
  }
  newHist.Rebin(2);
  return newHist;
}

//____________________________________________________________________72
const pair<vector<NuMatrixOutput*>,NuMatrixOutput*>
NuMatrixFitter::DoMultiRunPRLCCFit(const vector<NuMatrixInput*> mmInput,
                                   const vector<TH1D> numufdData,
                                   const NuModel_t model)
{
  vector<NuMatrixMethod*> numuExtrapolators;
  for (vector<NuMatrixInput*>::const_iterator inputIt = mmInput.begin();
       inputIt != mmInput.end();
       ++inputIt){
    numuExtrapolators.push_back
      (new NuMatrixMethod(**inputIt,
                          NuMMExtrapolation::kModularPRLCCFit));
  }
  
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  Double_t bestNorm = 0.0;

//   fsn2NuLow = 0.0;
//     fsn2NuHigh = 0.000;
//     fsn2NuGranularity = 0.001;
//     fdm2NuLow = 0.0e-3;
//     fdm2NuHigh = 0.0e-3;
//     fdm2NuGranularity = 0.005e-3;
    Double_t normLow =1.0;
    Double_t normHigh = 1.0;
    Double_t normGranularity = 0.0005;

  
  //Create the Chi2 histogram
  Int_t numSn2Bins =
    (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
  Int_t numDm2Bins =
    (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
  TH2D hChi2("hChi2",
             "",
             numSn2Bins,
             fsn2NuLow-fsn2NuGranularity/2.0,
             fsn2NuHigh+fsn2NuGranularity/2.0,
             numDm2Bins,
             fdm2NuLow-fdm2NuGranularity/2.0,
             fdm2NuHigh+fdm2NuGranularity/2.0); 
  
  for (Double_t sn2 = fsn2NuLow;
       sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
       sn2 += fsn2NuGranularity){
    for (Double_t dm2 = fdm2NuLow;
         dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
         dm2 += fdm2NuGranularity){
      
      Double_t chi2BestNorm = 1.0e10;
      Double_t bestNormThisLoop = 0.0;

      for (Double_t normalisation = normLow;
           normalisation <= normHigh+normGranularity/2.0;
           normalisation += normGranularity){
        
        Double_t chi2 = 0.0;
        
        vector<TH1D>::const_iterator fdDataIt = numufdData.begin();
        vector<NuMatrixMethod*>::iterator mmNuMuIt = numuExtrapolators.begin();
        for (;
             mmNuMuIt != numuExtrapolators.end();
             ++mmNuMuIt, ++fdDataIt){
          const TH1D* numufdPred = 0;
          if (kOscillation == model){
            numufdPred =
              (*mmNuMuIt)->PredictFDSpectrumBackgroundNoOsc(dm2,
                                                            sn2);
          }
          else if (kDecay == model){
            numufdPred =
              (*mmNuMuIt)->PredictFDDecaySpectrumBackgroundNoOsc(dm2,
                                                                 sn2);
          }
          else if (kDecoherence == model){
            numufdPred =
              (*mmNuMuIt)->PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
                                                                       sn2);
          }
          else {
            MSG("NuMatrixFitter.cxx",Msg::kError)
              << "Invalid neutrino disappearance model" << endl;
          }


          TH1D numufdPredScaled = this->RebinForFit(numufdPred);
          numufdPredScaled.Scale(normalisation);
          //      TH1D* pPredScaled = numufdPredScaled.Rebin(numFitBins,"ScaledPred",fitBinsArray);
          TH1D rebinnedData = this->RebinForFit(&(*fdDataIt));
//        TH1* pRebinnedData = rebinnedData.Rebin(numFitBins,"ScaledData",fitBinsArray);
          chi2 += this->CalculateLikelihood(&numufdPredScaled,&rebinnedData);
        }
        chi2 += this->NormalisationPenaltyTerm(normalisation);
        if (chi2 < bestChi2){
          bestChi2 = chi2;
          bestsn2 = sn2;
          bestdm2 = dm2;
          bestNorm = normalisation;
        }//if
        if (chi2 < chi2BestNorm){
          chi2BestNorm = chi2;
          bestNormThisLoop = normalisation;
        }//if
        MSG("NuMatrixFitter.cxx",Msg::kInfo)
          << "sn2 = " << sn2
          << ", dm2 = " << dm2
          << ", norm = " << normalisation
          << ", penalty = " << this->NormalisationPenaltyTerm(normalisation)
          << "; bestChi2 " << bestChi2 
          << "; bestsn2: " << bestsn2
          << "; bestdm2: " << bestdm2 
          << "; bestNorm: " << bestNorm << endl;
      }//normalisation
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      hChi2.SetBinContent(bin2D,chi2BestNorm);
    }//dm2
  }//sn2

  MSG("NuMatrixFitter.cxx",Msg::kInfo)
                << "Best fit:" << endl
                << "    sin2 = " << bestsn2 << endl
                << "    dm2 = " << bestdm2 << endl
                << "    norm = " << bestNorm << endl
                << "    chi2 = " << bestChi2 << endl;
  
  TH1D* hbestFitPred = 0;
  TH1D* hnoOscPred = 0;
  TH1D* hcombinedFDData = 0;
  vector<NuMatrixOutput*> mmOuts;
  Bool_t firstTime = true;
  vector<TH1D>::const_iterator fdDataIt = numufdData.begin();
  vector<NuMatrixMethod*>::iterator mmNuMuIt = numuExtrapolators.begin();
  for (;
       mmNuMuIt != numuExtrapolators.end();
       ++mmNuMuIt, ++fdDataIt){
    NuMatrixOutput* mmOutput = new NuMatrixOutput();
    mmOutput->NuMuChi2Surface(hChi2);
    mmOutput->NuMuFDData(*fdDataIt);
    mmOutput->NuMuBestFitFDPrediction(*(*mmNuMuIt)->
                                      PredictFDSpectrumBackgroundNoOsc(bestdm2,
                                                                       bestsn2));
    mmOutput->NuMuNoOscPrediction(*(*mmNuMuIt)->
                                  PredictFDSpectrumBackgroundNoOsc(0.0,
                                                                   0.0));
    mmOutput->BestFitPoint(bestsn2, bestdm2);
    mmOuts.push_back(mmOutput);
    if (firstTime){
      if (kOscillation == model){
        hbestFitPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDSpectrumBackgroundNoOsc(bestdm2,
                                            bestsn2));
        hnoOscPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDSpectrumBackgroundNoOsc(0.0,
                                            0.0));
      }
      else if (kDecay == model){
        hbestFitPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
                                                 bestsn2));
        hnoOscPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDDecaySpectrumBackgroundNoOsc(0.0,
                                                 0.0));
      }
      else if (kDecoherence == model){
        hbestFitPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
                                                 bestsn2));
        hnoOscPred = new TH1D
          (*(*mmNuMuIt)->
           PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
                                                 0.0));
      }
      hcombinedFDData = new TH1D(*fdDataIt);
      firstTime = false;
    }
    else{
      if (kOscillation == model){
        hbestFitPred->Add
          ((*mmNuMuIt)->
           PredictFDSpectrumBackgroundNoOsc(bestdm2,
                                            bestsn2));
        hnoOscPred->Add
          ((*mmNuMuIt)->
           PredictFDSpectrumBackgroundNoOsc(0.0,
                                            0.0));
      }
      else if (kDecay == model){
        hbestFitPred->Add
          ((*mmNuMuIt)->
           PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
                                                 bestsn2));
        hnoOscPred->Add
          ((*mmNuMuIt)->
           PredictFDDecaySpectrumBackgroundNoOsc(0.0,
                                                 0.0));
      }
      else if (kDecoherence == model){
        hbestFitPred->Add
          ((*mmNuMuIt)->
           PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
                                                 bestsn2));
        hnoOscPred->Add
          ((*mmNuMuIt)->
           PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
                                                 0.0));
      }
      hcombinedFDData->Add(&(*fdDataIt));
    }
  }
  
  NuMatrixOutput *mmOutputCombined = new NuMatrixOutput();
  mmOutputCombined->NuMuChi2Surface(hChi2);
  hcombinedFDData = new TH1D(this->RebinForFit(hcombinedFDData));
  mmOutputCombined->NuMuFDData(*hcombinedFDData);
  hbestFitPred = new TH1D(this->RebinForFit(hbestFitPred));
  hbestFitPred->Scale(bestNorm);
  mmOutputCombined->NuMuBestFitFDPrediction(*hbestFitPred);
  hnoOscPred = new TH1D(this->RebinForFit(hnoOscPred));
  mmOutputCombined->NuMuNoOscPrediction(*hnoOscPred);
  mmOutputCombined->BestFitPoint(bestsn2, bestdm2);

  pair<vector<NuMatrixOutput*>,NuMatrixOutput*> allMMOuts(mmOuts,mmOutputCombined);
  return allMMOuts;
}

//____________________________________________________________________72
void NuMatrixFitter::CPTFit(const string inputFilename,
                            const string fdDataFilename,
                            const Double_t dm2,
                            const Double_t sn2,
                            const string outputFilename)
{
  TFile fdDataFile(fdDataFilename.c_str(),"READ");
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  fdDataFile.Close();

  TFile *infile = new TFile(inputFilename.c_str(), "READ");
        

  NuMatrixInput *mmInput = 0;
  mmInput = (NuMatrixInput*) infile->Get("NuMatrixInput");
  if (!mmInput) {
    mmInput = new NuMatrixInput(infile);
  }
  

  
  TString objName = "NuMatrixOutput";
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
  }
  
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
    MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;

    //objName = fxmlConfig->FullName();
  }
  else {
    MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
    //objName = "Nominal";
  }

  NuMatrixOutput *mmOutput = this->DoCPTFit(mmInput, numubarfdData, dm2, sn2);
  if (fxmlConfig) {
    mmOutput->XMLConfig(fxmlConfig);
  }
  
  MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
  TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
  outputFile->cd();
  mmOutput->Write(objName);
  outputFile->Close();
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
}

//____________________________________________________________________72
void NuMatrixFitter::CCFitChargeCut(const string inputFilename,
                                    const string fdDataFilename,
                                    const string outputFilename,
                                    const Double_t nubarEcut,
                                    const NuModel_t model)
{  
  TFile fdDataFile(fdDataFilename.c_str(),"READ");
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  fdDataFile.Close();
  
  TFile *infile = new TFile(inputFilename.c_str(), "READ");
  
  NuMatrixInput mmInput(infile);
  
  TString objName = "NuMatrixOutput";
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
  }
  
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
    MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;

    //objName = xmlConfig->FullName();
  }
  else {
    MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
    //objName = "Nominal";
  }

  NuMatrixOutput *mmOutput = this->DoCCFitChargeCut(mmInput, 
                                                    numufdData, 
                                                    numubarfdData, 
                                                    nubarEcut,
                                                    model);

  if (fxmlConfig) {
    mmOutput->XMLConfig(fxmlConfig);
  }
  MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
  TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
  outputFile->cd();
  mmOutput->Write(objName);
  outputFile->Close();
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
}

//____________________________________________________________________72
void NuMatrixFitter::NoChargeCutFit(const string inputFilename,
                                    const string fdDataFilename,
                                    const string outputFilename,
                                    const NuModel_t model)
{  
  TFile fdDataFile(fdDataFilename.c_str(),"READ");
  TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyAll_FD");
  fdData->SetDirectory(0);
  fdData->Sumw2();
  fdDataFile.Close();
  
  TFile *infile = new TFile(inputFilename.c_str(), "READ");
  
  NuMatrixInput mmInput(infile);

  
  TString objName = "NuMatrixOutput";
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
  }
  
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
    MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;

    //objName = fxmlConfig->FullName();
  }
  else {
    MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
    //objName = "Nominal";
  }

  NuMatrixOutput *mmOutput = this->DoNoChargeCutFit(mmInput, fdData, model);

  if (fxmlConfig) {
    mmOutput->XMLConfig(fxmlConfig);
  }
  MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
  TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
  outputFile->cd();
  mmOutput->Write(objName);
  outputFile->Close();
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
}

//____________________________________________________________________72
void NuMatrixFitter::PRLCCFit(const string inputFilename,
                              const string fdDataFilename,
                              const string outputFilename,
                              const NuModel_t model)
{  
  MSG("NuMatrixFitter",Msg::kInfo)  << "NuMatrixFitter::PRLCCFit()" << endl;
  TFile fdDataFile(fdDataFilename.c_str(),"READ");
  TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  fdData->SetDirectory(0);
  fdData->Sumw2();
  fdDataFile.Close();
  
  TFile *infile = new TFile(inputFilename.c_str(), "READ");
  
  NuMatrixInput mmInput(infile);

  TString objName = "NuMatrixOutput";
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
  }
  
  MSG("NuMatrixFitter",Msg::kInfo)  << "fxmlConfig=" << fxmlConfig << endl;
  if (fxmlConfig) {
    FillGridSearchParams(fxmlConfig);
    MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;

    //objName = fxmlConfig->FullName();
  }
  else {
    MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
    //objName = "Nominal";
  }

  NuMatrixOutput *mmOutput = this->DoPRLCCFit(mmInput, fdData, model);
  if (fxmlConfig) {
     mmOutput->XMLConfig(fxmlConfig);
  }
  MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
  TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
  outputFile->cd();
  mmOutput->Write(objName);
  outputFile->Close();
  infile->cd();
  MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
  //    infile->Close();
}

//____________________________________________________________________72
void NuMatrixFitter::TransitionFit(const string inputFilename,
                                   const string fdDataFilename,
                                   const string outputFilename,
                                   const Int_t experiments)
{  
    TFile fdDataFile(fdDataFilename.c_str(),"READ");
    TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
    fdData->SetDirectory(0);
    fdData->Sumw2();
    fdDataFile.Close();
    
    TFile *infile = new TFile(inputFilename.c_str());
    TString objName = "NuMatrixOutput";
    if (!fxmlConfig) {
        fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
    }
    infile->Close();
    
    
    MSG("NuMatrixFitter",Msg::kInfo) << "fxmlConfig=" << fxmlConfig << endl;
    if (fxmlConfig) {
        FillGridSearchParams(fxmlConfig);
        MSG("NuMatrixFitter",Msg::kInfo) << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
        
        //objName = fxmlConfig->FullName();
    }
    else {
        MSG("NuMatrixFitter",Msg::kInfo) << "No NuXMLConfig object -- nominal fit assumed." << endl;
        //objName = "Nominal";
    }
    
    
    MSG("NuMatrixFitter",Msg::kInfo) << "Recreating file " << outputFilename << endl;
    TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");

    NuFluctuator *fl = new NuFluctuator();
    TH1D *fd_sub_Data = 0;
    TH1F *bestFits = new TH1F("bestFits","Best Fit Points",40,0.0,1.0);
    NuMatrixOutput *mmOutput = DoTransitionFit(inputFilename.c_str()/*mmInput*/, 
                                               fdData, 
                                               fxmlConfig->DM2Nu(), 
                                               fxmlConfig->SN2Nu());
    
    if (fxmlConfig) mmOutput->XMLConfig(fxmlConfig);
    
    TString name = objName;
    name += "_";
    name += "nofluct";
    outputFile->WriteObject(mmOutput,name);
    delete fd_sub_Data;
    
    for (int i = 0; i < experiments; i++) {
        fd_sub_Data = new TH1D(fl->Fluctuate(fdData));
        NuMatrixOutput *mmOutput = DoTransitionFit(inputFilename.c_str()/*mmInput*/, 
                                                   fd_sub_Data, 
                                                   fxmlConfig->DM2Nu(), 
                                                   fxmlConfig->SN2Nu());
        
        if (fxmlConfig) mmOutput->XMLConfig(fxmlConfig);

        bestFits->Fill(mmOutput->bestTransitionProb);
        TString name = objName;
        name += "_";
        name += i;
        outputFile->WriteObject(mmOutput,name);
        delete fd_sub_Data;
    }
    outputFile->WriteObject(bestFits,bestFits->GetName());
    outputFile->Close();
    MSG("NuMatrixFitter",Msg::kInfo) << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
}

//____________________________________________________________________72
void NuMatrixFitter::AppearanceAnalysis()
{
  //Get the FD data & PoT histograms
  //Data:
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  //PoT:
  TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
  hfdPoTTorTGT->SetDirectory(0);
  TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
  hfdSpillsPerFile->SetDirectory(0);
  TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
  hfdRun->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
  hndPoTTorTGT->SetDirectory(0);
  ndDataFile.Close();
  
  //Calculate the total PoT
  Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
  Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
  //   Double_t fdPoT = 6.5e20;
  MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;
  
  //Hard coded, unchecked fiducial masses. Should make configurable.
  Double_t ndFidMass = 45.128*Munits::tonne;
  Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
  
  NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMu.root");
  mmNuMu.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
  
  NuMatrixMethod mmNuMuBar(NuBinningScheme::kUnknown,
                           ndPoT,
                           fdPoT,
                           ffluxPoT,
                           ndFidMass,
                           fdFidMass,
                           fhelperFilename,
                           fxsecFilename,
                           fndDataFilename,
                           "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMuBar.root");
  mmNuMuBar.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
  
  mmNuMu.PredictFDFluxUnosc();
  mmNuMuBar.PredictFDFluxUnosc();
  
  mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
  mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());
  Double_t dm2 = 2.5e-3;
  Double_t sn2 = 1.0;
  Double_t dm2bar = 0.0;//2.5e-3;
  Double_t sn2bar = 0.0;//1.0;
  mmNuMu.SetFDAppearedFidFlux(mmNuMuBar.GetFDPotentialAppearanceFidFlux());
  mmNuMuBar.SetFDAppearedFidFlux(mmNuMu.GetFDPotentialAppearanceFidFlux());
  
  
  //  const TH1D* numufdPred =
  mmNuMu.PredictFDSpectrumAppearance(dm2,
                                     sn2,
                                     dm2bar,
                                     sn2bar,
                                     0.2);
  //  const TH1D* numubarfdPred =
  mmNuMuBar.PredictFDSpectrumAppearance(dm2bar,
                                        sn2bar,
                                        dm2,
                                        sn2,
                                        0.2);
  
  
  TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMu.root","RECREATE");
  numufdData->Write();
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}
  
  TFile *sfilebar = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMuBar.root","RECREATE");
  numubarfdData->Write();
  sfilebar->Close();
  if (sfilebar){delete sfilebar; sfilebar = 0;}
}

//____________________________________________________________________72
void NuMatrixFitter::CPTAnalysis()
{
  //Get the FD data & PoT histograms
  //Data:
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  //PoT:
  TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
  hfdPoTTorTGT->SetDirectory(0);
  TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
  hfdSpillsPerFile->SetDirectory(0);
  TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
  hfdRun->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
  hndPoTTorTGT->SetDirectory(0);
  ndDataFile.Close();

  //Calculate the total PoT
  Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
   Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
   //   Double_t fdPoT = 6.5e20;
  MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;

  //Hard coded, unchecked fiducial masses. Should make configurable.
  Double_t ndFidMass = 45.128*Munits::tonne;
  Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);

  NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");
  mmNuMu.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);

  NuMatrixMethod mmNuMuBar(NuBinningScheme::kUnknown,
                           ndPoT,
                           fdPoT,
                           ffluxPoT,
                           ndFidMass,
                           fdFidMass,
                           fhelperFilename,
                           fxsecFilename,
                           fndDataFilename,
                           "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMuBar.root");
  mmNuMuBar.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);

  mmNuMu.PredictFDFluxUnosc();
  mmNuMuBar.PredictFDFluxUnosc();

  mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
  mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());

  // TODO: When we start using this code again, use the member
  // variables for the grid search parameters
  Double_t bestChi2 = 1.0e10;
  Double_t bestsn2 = 0.0;
  Double_t bestdm2 = 0.0;
  Double_t bestsn2bar = 0.0;
  Double_t bestdm2bar = 0.0;
  Double_t sn2 = 0;
  Double_t dm2 = 0;
  Double_t sn2bar = 0;
  Double_t dm2bar = 0;
  Double_t sn2Low = 0.8;
  Double_t sn2High = 1.0;
  Double_t sn2Granularity = 0.005;
  Double_t dm2Low = 1.8e-3;
  Double_t dm2High = 3.2e-3;
  Double_t dm2Granularity = 0.05e-3;
  Double_t sn2BarLow = 0.0;
  Double_t sn2BarHigh = 1.0;
  Double_t sn2BarGranularity = 0.01;
  Double_t dm2BarLow = 0e-3;
  Double_t dm2BarHigh = 10e-3;
  Double_t dm2BarGranularity = 0.1e-3;
  /*
  for (sn2 = sn2Low; sn2 <= sn2High; sn2 += sn2Granularity){
    for (dm2 = dm2Low; dm2 <= dm2High; dm2 += dm2Granularity){
      for (sn2bar = sn2BarLow; sn2bar <= sn2BarHigh; sn2bar += sn2BarGranularity){
        for (dm2bar = dm2BarLow; dm2bar <= dm2BarHigh; dm2bar += dm2BarGranularity){
          const TH1D* numufdPred =
            mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                          sn2,
                                                          dm2bar,
                                                          sn2bar);
          const TH1D* numubarfdPred =
            mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                             sn2bar,
                                                             dm2,
                                                             sn2);
          Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
          chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
          MSG("NuMatrixFitter",Msg::kInfo)  << "dm2 " << dm2
               << "; sn2 " << sn2
               << "; dm2bar " << dm2bar
               << "; sn2bar " << sn2bar
               << "; chi2 " << chi2
               << endl;
          if (chi2 < bestChi2){
            bestChi2 = chi2;
            bestsn2 = sn2;
            bestdm2 = dm2;
            bestsn2bar = sn2bar;
            bestdm2bar = dm2bar;
          }//if
        }//dm2bar
      }//sn2bar
    }//dm2
  }//sn2
  */
 
  MSG("NuMatrixFitter",Msg::kInfo)  << "Best fit: dm2 = " << bestdm2
       << "; sn2 = " << bestsn2
       << "; dm2bar = " << bestdm2bar
       << "; sn2bar = " << bestsn2bar
       << "; bestChi2 = " << bestChi2 << endl;

  Int_t numSn2bins = (Int_t) round((sn2High-sn2Low)/sn2Granularity) + 1;
  Int_t numDm2bins = (Int_t) round((dm2High-dm2Low)/dm2Granularity) + 1;
  //Loop again to fill chi2 surfaces.
  TH2D hChi2("hChi2",
             "",
             numSn2bins,
             sn2Low-sn2Granularity/2.0,
             sn2High+sn2Granularity/2.0,
             numDm2bins,
             dm2Low-dm2Granularity/2.0,
             dm2High+dm2Granularity/2.0);
  
//   sn2bar = bestsn2bar;
//   dm2bar = bestdm2bar;
  sn2bar = 1.0;
  dm2bar = 6.0e-3;
  
  for (sn2 = sn2Low;
       sn2 < sn2High+sn2Granularity/2.0;
       sn2 += sn2Granularity){
    for (dm2 = dm2Low;
         dm2 < dm2High+dm2Granularity/2.0;
         dm2 += dm2Granularity){
      const TH1D* numufdPred =
        mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                      sn2,
                                                      dm2bar,
                                                      sn2bar);
      const TH1D* numubarfdPred =
        mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                         sn2bar,
                                                         dm2,
                                                         sn2);
      Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
      chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
      Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
      Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
      Int_t bin2D = hChi2.GetBin(xBin,yBin);
      MSG("NuMatrixFitter",Msg::kInfo)  << "Setting bin " << sn2
           << " and " << dm2
           << " with " << chi2
           << endl;
      hChi2.SetBinContent(bin2D,chi2);
    }//dm2
  }//sn2
  
  Int_t numSn2BarBins =
    (Int_t) round((sn2BarHigh-sn2BarLow)/sn2BarGranularity) + 1;
  Int_t numDm2BarBins =
    (Int_t) round((dm2BarHigh-dm2BarLow)/dm2BarGranularity) + 1;
  //Loop again to fill chi2bar surfaces.
  
  TH2D hChi2Bar("hChi2Bar",
             "",
             numSn2BarBins,
             sn2BarLow-sn2BarGranularity/2.0,
             sn2BarHigh+sn2BarGranularity/2.0,
             numDm2BarBins,
             dm2BarLow-dm2BarGranularity/2.0,
             dm2BarHigh+dm2BarGranularity/2.0);
  
//   sn2 = bestsn2;
//   dm2 = bestdm2;
  sn2 = 1.0;
  dm2 = 2.3e-3;
  for (sn2bar = sn2BarLow;
       sn2bar < sn2BarHigh+sn2BarGranularity/2.0;
       sn2bar += sn2BarGranularity){
    for (dm2bar = dm2BarLow;
         dm2bar < dm2BarHigh+dm2BarGranularity/2.0;
         dm2bar += dm2BarGranularity){
      const TH1D* numufdPred =
        mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                      sn2,
                                                      dm2bar,
                                                      sn2bar);
      const TH1D* numubarfdPred =
        mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                         sn2bar,
                                                         dm2,
                                                         sn2);
      Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
      chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
      Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
      Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
      MSG("NuMatrixFitter",Msg::kInfo)  << "Setting bin " << sn2bar
           << " and " << dm2bar
           << " with " << chi2
           << endl;
      Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
      hChi2Bar.SetBinContent(bin2D,chi2);
    }//dm2Bar
  }//sn2Bar
  
  
  TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root","RECREATE");
  numufdData->Write();
  hChi2.Write();
  hChi2Bar.Write();
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}

  TFile *sfilebar = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMuBar.root","RECREATE");
  numubarfdData->Write();
  sfilebar->Close();
  if (sfilebar){delete sfilebar; sfilebar = 0;}
}

//____________________________________________________________________72
void NuMatrixFitter::WriteHistosForFitter(TString outputFileName)
{

  //Get the FD data & PoT histograms
  //Data:
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  //PoT:
  TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
  hfdTotalPot->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
  hndTotalPot->SetDirectory(0);
  
  //Open the helper file.
  TFile helperFile(fhelperFilename.c_str(),"READ");

  // Get the NuXMLConfig Object. Look in the ND data file first. If
  // it's not in there, look in the helper file.
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
  }
  if (!fxmlConfig){
     fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
   }

  ndDataFile.Close();

  //Calculate the total PoT
  Double_t ndPoT = hndTotalPot->Integral();
  Double_t fdPoT = hfdTotalPot->Integral();
  if (fxmlConfig->OverriddenFDPoTs()>0) {
    fdPoT = fxmlConfig->OverriddenFDPoTs();
  }

  // If scaling is used, the scaling should have been done by this
  // stage, so use the scaled pots instead of the overridden pots
  if (fxmlConfig->ScaledFDPoTs()>0) {
    fdPoT = fxmlConfig->ScaledFDPoTs();
  }

  MSG("NuMatrixFitter",Msg::kInfo) << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;

  NuCuts nuCuts;
  Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
                                           SimFlag::kMC,
                                           NuCuts::kJJE1);
  Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
                                           SimFlag::kMC,
                                           NuCuts::kJJE1);
  MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
       << "; far mas: " << fdFidMass << endl;

  //Get the binning scheme to use
  NuBinningScheme::NuBinningScheme_t binningScheme =
    static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
  MSG("NuMatrixFitter",Msg::kInfo)  << "Using binningScheme " << binningScheme << endl;

  NuMatrixMethod mmNuMu(binningScheme,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        "");

  //pull in all the histograms from the file
  mmNuMu.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);

  //set all the parameters
  //pass in the filenames
  NuMatrixMethod mmNuMuBar(binningScheme,
                           ndPoT,
                           fdPoT,
                           ffluxPoT,
                           ndFidMass,
                           fdFidMass,
                           fhelperFilename,
                           fxsecFilename,
                           fndDataFilename,
                           "");

  //pull in all the histograms from the file
  mmNuMuBar.SetExtrapolationScheme
    (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
  cout << "Debug Start" << endl;
  //do lots and lots of steps
  mmNuMu.PredictFDFluxUnosc();
  mmNuMuBar.PredictFDFluxUnosc();
  cout << "Debug End" << endl;
  //use extrapolated numu spectrum to predict numu contamination in 
  //numubar spectrum (and vice versa)
  //The "K" spectrum to pass out to external fitter is now created
  mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
  mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());

  //For a nubar appearance analyses: Use the exptrapolated numu
  //spectrum to create the spectrum of nubars that could appear. This
  //then needs to be oscillated (numu DISAPPEARANCE) to get the actual
  //spectrum that could appear. Oscillation is done in the fitting code.
  mmNuMu.SetFDAppearedFidFlux(mmNuMuBar.GetFDPotentialAppearanceFidFlux());
  mmNuMuBar.SetFDAppearedFidFlux(mmNuMu.GetFDPotentialAppearanceFidFlux());

  mmNuMu.PredictFDSpectrumBackgroundOscCombined(0.0025,1,
                                                0.0025,1);
  mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(0.0025,1,
                                                   0.0025,1);


  //write out files
  //create the file
  TFile* sfile = new TFile(outputFileName,"RECREATE");
  MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
  sfile->Print();

  fxmlConfig->Write();
  
  //write files
  mmNuMu.WriteFilesForFitter("NM");
  mmNuMuBar.WriteFilesForFitter("NMB");

  //close file
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}  


/*
          const TH1D* numufdPred =
            mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                          sn2,
                                                          dm2bar,
                                                          sn2bar);
          const TH1D* numubarfdPred =
            mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                             sn2bar,
                                                             dm2,
                                                             sn2);
*/

}

//____________________________________________________________________72
void NuMatrixFitter
::WriteMultiRunHistosForFitter(const vector<TString> outputFileNames)
{
  //Hard coded, unchecked fiducial masses. Should make configurable.
  Double_t ndFidMass = 45.128*Munits::tonne;
  Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()/52.92;

  //FD data
  vector<TH1D> numufdData;
  vector<TH1D> nubarfdData;
  //NuMu extrapolators
  vector<NuMatrixMethod*> numuExtrapolators;
  vector<NuMatrixMethod*> nubarExtrapolators;
  //XML configs
  vector<NuXMLConfig> xmlConfigs;
  
  vector<string>::const_iterator helperFileIt = fvhelperFilenames.begin();
  vector<string>::const_iterator fdFileIt = fvfdDataFilenames.begin();
  for (vector<string>:: const_iterator ndFileIt = fvndDataFilenames.begin();
       ndFileIt != fvndDataFilenames.end();
       ++ndFileIt, ++fdFileIt, ++helperFileIt){
    //Calculate ND PoT
    TFile ndDataFile((*ndFileIt).c_str(),"READ");
    TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
    Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
    ndDataFile.Close();

    //Get FD data
    TFile fdDataFile((*fdFileIt).c_str(),"READ");
    TH1D* numuHist = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
    numuHist->SetDirectory(0);
    numuHist->Sumw2();
    numufdData.push_back(*numuHist);
    TH1D* numubarHist = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
    numubarHist->SetDirectory(0);
    numubarHist->Sumw2();
    nubarfdData.push_back(*numubarHist);

    //Get FD PoT
    TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
    TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
    Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
    fdDataFile.Close();

    //Get NuXMLConfig objects
    TFile helperFile((*helperFileIt).c_str(),"READ");
    NuXMLConfig *xmlConfig = (NuXMLConfig*) helperFile.Get("NuXMLConfig");
    xmlConfigs.push_back(*xmlConfig);

    //Binning scheme
    NuBinningScheme::NuBinningScheme_t binningScheme =
      static_cast<NuBinningScheme::NuBinningScheme_t>(xmlConfig->BinningScheme());
    MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme" << binningScheme << endl;

    //Construct the extrapolators
    NuMatrixMethod* numuExtrap = new NuMatrixMethod(binningScheme,
                                                    ndPoT,
                                                    fdPoT,
                                                    ffluxPoT,
                                                    ndFidMass,
                                                    fdFidMass,
                                                    *helperFileIt,
                                                    fxsecFilename,
                                                    *ndFileIt,
                                                    "");
    numuExtrap->SetExtrapolationScheme
      (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
    numuExtrapolators.push_back(numuExtrap);
    
    NuMatrixMethod* nubarExtrap = new NuMatrixMethod(binningScheme,
                                                     ndPoT,
                                                     fdPoT,
                                                     ffluxPoT,
                                                     ndFidMass,
                                                     fdFidMass,
                                                     *helperFileIt,
                                                     fxsecFilename,
                                                     *ndFileIt,
                                                     "");
    nubarExtrap->SetExtrapolationScheme
      (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
    nubarExtrapolators.push_back(nubarExtrap);
  }

  //Now do the extrapolation steps
  vector<NuMatrixMethod*>::iterator numuExtrapIt = numuExtrapolators.begin();
  vector<NuMatrixMethod*>::iterator nubarExtrapIt = nubarExtrapolators.begin();
  vector<TString>::const_iterator outfileIt = outputFileNames.begin();
  vector<NuXMLConfig>::iterator xmlConfigIt = xmlConfigs.begin();
  for (;
       numuExtrapIt != numuExtrapolators.end();
       ++numuExtrapIt, ++nubarExtrapIt, ++outfileIt){
    (*numuExtrapIt)->PredictFDFluxUnosc();
    (*nubarExtrapIt)->PredictFDFluxUnosc();

    (*numuExtrapIt)->SetFDCCTrueUnoscBackground
      ((*nubarExtrapIt)->GetFDOtherCCTrueUnoscBackground());
    (*nubarExtrapIt)->SetFDCCTrueUnoscBackground
      ((*numuExtrapIt)->GetFDOtherCCTrueUnoscBackground());

    //For nubar appearance analysis
    (*numuExtrapIt)->SetFDAppearedFidFlux
      ((*nubarExtrapIt)->GetFDPotentialAppearanceFidFlux());
    (*nubarExtrapIt)->SetFDAppearedFidFlux
      ((*numuExtrapIt)->GetFDPotentialAppearanceFidFlux());

    //Don't think we need to run this stage:
    (*numuExtrapIt)->PredictFDSpectrumBackgroundOscCombined(0.0025,1,
                                                          0.0025,1);
    (*nubarExtrapIt)->PredictFDSpectrumBackgroundOscCombined(0.0025,1,
                                                           0.0025,1);

    //Write the output files
    TFile* sfile = new TFile(*outfileIt,"RECREATE");
    MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
    sfile->Print();
    (*xmlConfigIt).Write();
    (*numuExtrapIt)->WriteFilesForFitter("NM");
    (*nubarExtrapIt)->WriteFilesForFitter("NMB");
    sfile->Close();
    if (sfile){delete sfile; sfile = 0;} 
  }
}

//____________________________________________________________________72

void NuMatrixFitter
::WriteNoChargeCutHistosForFitter(std::string outputFileName)
{
  //Get the FD data & PoT histograms
  //Data:
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  //PoT:
  TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
  hfdTotalPot->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
  hndTotalPot->SetDirectory(0);
  
  //Open the helper file.
  TFile helperFile(fhelperFilename.c_str(),"READ");

  // Get the NuXMLConfig Object. Look in the ND data file first. If
  // it's not in there, look in the helper file.
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
  }
  if (!fxmlConfig){
    fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
  }

  ndDataFile.Close();

  //Calculate the total PoT
  Double_t ndPoT = hndTotalPot->Integral();
  Double_t fdPoT = hfdTotalPot->Integral();
  if (fxmlConfig->OverriddenFDPoTs()>0) {
    fdPoT = fxmlConfig->OverriddenFDPoTs();
  }

  // If scaling is used, the scaling should have been done by this
  // stage, so use the scaled pots instead of the overridden pots
  if (fxmlConfig->ScaledFDPoTs()>0) {
    fdPoT = fxmlConfig->ScaledFDPoTs();
  }

  MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;

  NuCuts nuCuts;
  Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
                                           SimFlag::kMC,
                                           NuCuts::kJJE1);
  Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
                                           SimFlag::kMC,
                                           NuCuts::kJJE1);

  MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
       << "; far mas: " << fdFidMass << endl;

  //Binning scheme
  NuBinningScheme::NuBinningScheme_t binningScheme =
    static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
  MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme " << binningScheme << endl;

  NuMatrixMethod mmNuMu(binningScheme,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");

  //pull in all the histograms from the file
  mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNoChargeCut);
  mmNuMu.PredictFDFluxUnosc();
  

  //write out files
  //create the file
  TFile* sfile=0;
  if (outputFileName=="") {
    sfile=new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMNoCutOutputForFitter.root","RECREATE");
  }
  else {
    sfile=new TFile(outputFileName.c_str(),"RECREATE");
  }
  MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
  sfile->Print();
  
  fxmlConfig->Write();
  
  //write files
  mmNuMu.WriteFilesForFitter("All");
  
  //close file
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}  


/*
          const TH1D* numufdPred =
            mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
                                                          sn2,
                                                          dm2bar,
                                                          sn2bar);
          const TH1D* numubarfdPred =
            mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
                                                             sn2bar,
                                                             dm2,
                                                             sn2);
*/
}

//____________________________________________________________________72
void NuMatrixFitter
::WritePRLCCHistosForFitter(std::string outputFileName)
{
  //Get the FD data & PoT histograms
  //Data:
  TFile fdDataFile(ffdDataFilename.c_str(),"READ");
  TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
  numufdData->SetDirectory(0);
  numufdData->Sumw2();
  /*
  TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
  numubarfdData->SetDirectory(0);
  numubarfdData->Sumw2();
  */
  //PoT:
  TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
  hfdTotalPot->SetDirectory(0);
  fdDataFile.Close();

  //Get the ND PoT histograms
  TFile ndDataFile(fndDataFilename.c_str(),"READ");
  TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
  hndTotalPot->SetDirectory(0);
  
  //Open the helper file.
  TFile helperFile(fhelperFilename.c_str(),"READ");

  // Get the NuXMLConfig Object. Look in the ND data file first. If
  // it's not in there, look in the helper file.
  if (!fxmlConfig) {
    fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
  }
  if (!fxmlConfig){
    fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
  }

  ndDataFile.Close();

  //Calculate the total PoT
  Double_t ndPoT = hndTotalPot->Integral();
  Double_t fdPoT = hfdTotalPot->Integral();
  if (fxmlConfig->OverriddenFDPoTs()>0) {
    fdPoT = fxmlConfig->OverriddenFDPoTs();
  }

  // If scaling is used, the scaling should have been done by this
  // stage, so use the scaled pots instead of the overridden pots
  if (fxmlConfig->ScaledFDPoTs()>0) {
    fdPoT = fxmlConfig->ScaledFDPoTs();
  }

  MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
       << "; fdPoT: " << fdPoT
       << endl;

  NuCuts nuCuts;
  Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
                                           SimFlag::kData,
                                           NuCuts::kJJE1);
  Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
                                           SimFlag::kData,
                                           NuCuts::kJJE1);
  
  MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
       << "; far mas: " << fdFidMass << endl;

  //Binning scheme
  NuBinningScheme::NuBinningScheme_t binningScheme =
    static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
  MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme " << binningScheme << endl;

  NuMatrixMethod mmNuMu(binningScheme,
                        ndPoT,
                        fdPoT,
                        ffluxPoT,
                        ndFidMass,
                        fdFidMass,
                        fhelperFilename,
                        fxsecFilename,
                        fndDataFilename,
                        "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");

  //pull in all the histograms from the file
  mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kPRLCCAnalysis);
  mmNuMu.PredictFDFluxUnosc();
  

  //write out files
  //create the file
  TFile* sfile=0;
  if (outputFileName=="") {
    sfile=new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMNoCutOutputForFitter.root","RECREATE");
  }
  else {
    sfile=new TFile(outputFileName.c_str(),"RECREATE");
  }
  MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
  sfile->Print();
  
  fxmlConfig->Write();
  
  //write files
  mmNuMu.WriteFilesForFitter("NM");
  
  //close file
  sfile->Close();
  if (sfile){delete sfile; sfile = 0;}  

}

//____________________________________________________________________72

void NuMatrixFitter::NuMuBarAppearanceAnalysis()
{

}

---