NuXFitAnalysis.cxx

Go to the documentation of this file.

// $Id: NuXFitAnalysis.cxx,v 1.13 2008/05/02 19:39:22 evans Exp $
//
// Combined numu & numubar extrapolation and FD fit by cross-section
// fitting
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
//
// $Log: NuXFitAnalysis.cxx,v $
// Revision 1.13  2008/05/02 19:39:22  evans
// Updating a comment.
//
// Revision 1.12  2008/01/17 22:34:41  evans
// ...and now the transition analysis is configurable between Minuit and
// grid search (well, line search at least).
//
// Revision 1.11  2008/01/17 22:23:06  evans
// CPT fit is now configurable between grid search and Minuit.
//
// Revision 1.10  2008/01/11 00:12:52  hartnell
//
// Have to be careful not to confuse hist->Reset() with hist->Clear().
//
// Change comment from:
//
// //Clear recoE histograms
//
// to
//
// //Reset recoE histograms (removes bin contents and errors)
//
// to reflect what is done. Just a little paranoia...
//
// Revision 1.9  2008/01/10 17:58:33  evans
// Adding IsGoodMajorityCurvature to the list of selection cuts.
//
// Revision 1.8  2007/12/27 22:13:46  evans
// New class to allow the XFit extrapolation to be performed on multiple
// runs (e.g. SKZP periods), along with necessary interface additions to
// NuXFitAnalysis.
//
// Revision 1.7  2007/12/25 21:12:56  evans
// In the transition analysis, using the true neutrino energy for the
// efficiency corrections instead of the reconstructed energy.
//
// Revision 1.6  2007/12/21 17:41:18  evans
// Debugging and creating a FD output file for the XFit CPT analysis.
//
// Revision 1.5  2007/12/20 19:12:45  evans
// Enabling the XFit analysis to receive its input data as a
// reconstructed histogram (as well as a NuDST: you can do either). This
// means it can fit the output of the Fluctuator in either the ND or FD
// (or both).
//
// Revision 1.4  2007/12/07 14:53:15  evans
// Adding protection against uninitialised variables found by Robert
// Hatcher when compiling with gcc+optimization.
//
// Revision 1.3  2007/12/04 13:13:47  evans
// Adding a mechanism to return the best fit parameters.
//
// Revision 1.2  2007/12/03 17:10:19  evans
// Adding code to do a numu->numubar transition analysis withtruth
// fitting. Also added the ability to do a stystematics study with the
// NuSystematic class.
//
// Revision 1.1  2007/11/15 13:46:27  evans
// The class to apply my ND XSec weights to the far detector data and do
// an FD CPT fit (using Minuit2).
//
//

#include "Minuit2/ContoursError.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/FunctionMinimum.h"
#include "TFile.h"
#include "TFitterMinuit.h"
#include "TGraph.h"
#include "TGraph2DErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TMinuit.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuTreeWrapper.h"
#include "NtupleUtils/NuXFitAnalysis.h"
#include "NtupleUtils/NuXMLConfig.h"

ClassImp(NuXFitAnalysis)

CVSID("$Id: NuXFitAnalysis.cxx,v 1.13 2008/05/02 19:39:22 evans Exp $");

using namespace NuXSecFit;
using namespace ROOT::Minuit2;

using std::string;
using std::vector;

//____________________________________________________________________72
NuXFitAnalysis::NuXFitAnalysis()
  : fdoGridSearch(false),
    fwriteOutput(true),
    fUp(1.0),
    ffdDataPoT(0.0),
    ffdMCPoT(0.0),
    fdm2(0.0),
    fsn2(0.0),
    fBestDm2(-1.0),
    fBestSn2(-1.0),
    fBestDm2Bar(-1.0),
    fBestSn2Bar(-1.0)
{
  fanalysisSetting = NuXFit::kUnknown;
  fanaVersion = NuCuts::kJJE1;
  fnuSystematic = new NuSystematic();

  fSelEffFileName = "";
  fxsecfilename = "";

  fOutputFile = 0;

  ftFDDataTree = 0;
  ftFDMCTree = 0;

  fhNuBarCCSelectionEfficiency = 0;
  fhNuMuCCSelectionEfficiency = 0;
  fNuBarCCCrossSectionGraph = 0;
  fNuMuCCCrossSectionGraph = 0;
  fhNuBarCCCrossSections = 0;
  fhNuMuCCCrossSections = 0;

  fhFDNuMuCCDataFromFile = 0;
  fhFDNuMuBarCCDataFromFile = 0;

  fhFDDataRecoENuMuCC = 0;
  fhFDDataRecoENuMuBarCC = 0;
  fhFDMCRecoENuMuCC = 0;
  fhFDMCRecoENuMuBarCC = 0;

  fFitter = new TFitterMinuit();
  fFitter->CreateMinimizer();
  fFitter->SetMinuitFCN(this);
}

//____________________________________________________________________72
NuXFitAnalysis::NuXFitAnalysis(const NuXMLConfig& xmlConfig)
  : fdoGridSearch(false),
    fUp(1.0),
    ffdDataPoT(0.0),
    ffdMCPoT(0.0),
    fdm2(0.0),
    fsn2(0.0)
{
  fanalysisSetting = NuXFit::kUnknown;
  fanaVersion = NuCuts::kJJE1;
  fnuSystematic = new NuSystematic(xmlConfig);

  fSelEffFileName = "";
  fxsecfilename = "";

  fOutputFile = 0;

  ftFDDataTree = 0;
  ftFDMCTree = 0;

  fhNuBarCCSelectionEfficiency = 0;
  fhNuMuCCSelectionEfficiency = 0;
  fNuBarCCCrossSectionGraph = 0;
  fNuMuCCCrossSectionGraph = 0;
  fhNuBarCCCrossSections = 0;
  fhNuMuCCCrossSections = 0;

  fhFDNuMuCCDataFromFile = 0;
  fhFDNuMuBarCCDataFromFile = 0;

  fhFDDataRecoENuMuCC = 0;
  fhFDDataRecoENuMuBarCC = 0;
  fhFDMCRecoENuMuCC = 0;
  fhFDMCRecoENuMuBarCC = 0;

  fFitter = new TFitterMinuit();
  fFitter->CreateMinimizer();
  fFitter->SetMinuitFCN(this);
}

//____________________________________________________________________72

NuXFitAnalysis::~NuXFitAnalysis()
{
  if (fOutputFile){
    fOutputFile->Close();
    delete fOutputFile;
    fOutputFile = 0;
  }

  if (fnuSystematic){delete fnuSystematic; fnuSystematic = 0;}

  if (fhFDNuMuCCDataFromFile){
    delete fhFDNuMuCCDataFromFile; fhFDNuMuCCDataFromFile = 0;
  }
  if (fhFDNuMuBarCCDataFromFile){
    delete fhFDNuMuBarCCDataFromFile; fhFDNuMuBarCCDataFromFile = 0;
  }
}

//____________________________________________________________________72
double NuXFitAnalysis
::operator () (const vector<double>& pars) const
{
  static Int_t fitcounter = 0;
  if (!(fitcounter%50)){
    cout << "Minuit FD XFit call " << fitcounter << endl;
  }
  ++fitcounter;
  if (NuXFit::kCPT == fanalysisSetting){
    this->FillCPTFDPrediction(pars);
  }
  else if(NuXFit::kTransition == fanalysisSetting){
    this->FillTransitionFDPrediction(pars);
  }
  else{
    MSG("NuXFitAnalysis.cxx",Msg::kFatal)
      << "Incorrect analysis setting" << endl;
    return -1.0;
  }
  return this->CalculateLikelihood();
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis::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);}
  }
//   cout << "Like = " << like << endl;
  return like;
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis::CalculateLikelihood() const
{
  //Aim to minimise -2lnL. That's what I'm returning.

  Double_t like = 0;

  for (Int_t i=1; i<=fhFDMCRecoENuMuCC->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fhFDMCRecoENuMuCC->GetBinContent(i);
    Double_t dnu = fhFDDataRecoENuMuCC->GetBinContent(i);
    if (mnu<0.0001){mnu = 0.0001;}
    if (dnu){
      like += 2*(mnu - dnu + dnu*log(dnu/mnu));
    }
    else{like += 2*(mnu - dnu);}
  }

  for (Int_t i=1; i<=fhFDMCRecoENuMuBarCC->GetNbinsX(); ++i){
    Double_t mbar = fhFDMCRecoENuMuBarCC->GetBinContent(i);
    Double_t dbar = fhFDDataRecoENuMuBarCC->GetBinContent(i);
    if (mbar<0.0001){mbar = 0.0001;}
    if (dbar){
      like += 2*(mbar - dbar + dbar*log(dbar/mbar));
    }
    else{like += 2*(mbar - dbar);}
  }
//   cout << "Like = " << like << endl;
  return like;
}

//____________________________________________________________________72
void NuXFitAnalysis::OutputFileName(const string outFile)
{
  if (fOutputFile){
    fOutputFile->Close();
    delete fOutputFile;
    fOutputFile = 0;
  }
  fOutputFile = new TFile(outFile.c_str(),"RECREATE");
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CreateFDNuMuBarCCDataHisto()
{
  if (fhFDDataRecoENuMuBarCC) {
    delete fhFDDataRecoENuMuBarCC; fhFDDataRecoENuMuBarCC = 0;
  }

  if (fhFDNuMuCCDataFromFile && fhFDNuMuBarCCDataFromFile){
    fhFDDataRecoENuMuBarCC = new TH1D(*fhFDNuMuBarCCDataFromFile);
    fhFDDataRecoENuMuBarCC->Reset();
    return true;
  }

  Int_t size = fvFDRecoENuMuBarCCBins.size();
  if (!size){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD reconstructed energy binning not supplied."
      << endl;
    return false;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         fvFDRecoENuMuBarCCBins.begin();
       it != fvFDRecoENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhFDDataRecoENuMuBarCC = new TH1D("fhFDDataRecoENuMuBarCC",
                               "FD NuMuBar CC data",
                               size-1,
                               recoEBinEdges);

  if (!fhFDDataRecoENuMuBarCC) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CreateFDNuMuCCDataHisto()
{
  if (fhFDDataRecoENuMuCC) {
    delete fhFDDataRecoENuMuCC; fhFDDataRecoENuMuCC = 0;
  }

  if (fhFDNuMuCCDataFromFile && fhFDNuMuBarCCDataFromFile){
    fhFDDataRecoENuMuCC = new TH1D(*fhFDNuMuCCDataFromFile);
    fhFDDataRecoENuMuCC->Reset();
    return true;
  }

  Int_t size = fvFDRecoENuMuCCBins.size();
  if (!size){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD reconstructed energy binning not supplied."
      << endl;
    return false;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         fvFDRecoENuMuCCBins.begin();
       it != fvFDRecoENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhFDDataRecoENuMuCC = new TH1D("fhFDDataRecoENuMuCC",
                               "FD NuMu CC data",
                               size-1,
                               recoEBinEdges);

  if (!fhFDDataRecoENuMuCC) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CreateFDNuMuBarCCMCHisto()
{
  if (fhFDMCRecoENuMuBarCC) {
    delete fhFDMCRecoENuMuBarCC; fhFDMCRecoENuMuBarCC = 0;
  }

  if (fhFDNuMuCCDataFromFile && fhFDNuMuBarCCDataFromFile){
    fhFDMCRecoENuMuBarCC = new TH1D(*fhFDNuMuBarCCDataFromFile);
    fhFDMCRecoENuMuBarCC->Reset();
    return true;
  }

  Int_t size = fvFDRecoENuMuBarCCBins.size();
  if (!size){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD reconstructed energy binning not supplied."
      << endl;
    return false;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         fvFDRecoENuMuBarCCBins.begin();
       it != fvFDRecoENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhFDMCRecoENuMuBarCC = new TH1D("fhFDMCRecoENuMuBarCC",
                               "FD NuMuBar CC prediction",
                               size-1,
                               recoEBinEdges);

  if (!fhFDMCRecoENuMuBarCC) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CreateFDNuMuCCMCHisto()
{
  if (fhFDMCRecoENuMuCC) {
    delete fhFDMCRecoENuMuCC; fhFDMCRecoENuMuCC = 0;
  }

  if (fhFDNuMuCCDataFromFile && fhFDNuMuBarCCDataFromFile){
    fhFDMCRecoENuMuCC = new TH1D(*fhFDNuMuCCDataFromFile);
    fhFDMCRecoENuMuCC->Reset();
    return true;
  }

  Int_t size = fvFDRecoENuMuCCBins.size();
  if (!size){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD reconstructed energy binning not supplied."
      << endl;
    return false;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         fvFDRecoENuMuCCBins.begin();
       it != fvFDRecoENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhFDMCRecoENuMuCC = new TH1D("fhFDMCRecoENuMuCC",
                               "FD NuMu CC prediction",
                               size-1,
                               recoEBinEdges);

  if (!fhFDMCRecoENuMuCC) return false;
  return true;
}

//____________________________________________________________________72
void NuXFitAnalysis::FillCPTFDPrediction
(const vector<double>& par) const
{
  //Reset recoE histograms (removes bin contents and errors)
  fhFDMCRecoENuMuCC->Reset();
  fhFDMCRecoENuMuBarCC->Reset();

  Double_t shwEShift=1.0;
  Double_t trkEShift = 1.0;

  Double_t Dm2 = par[0];
  Double_t Sn2 = par[1];
  Double_t Dm2Bar = par[2];
  Double_t Sn2Bar = par[3];

  //Fill the nubar histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itBar
         = fvFDMCCCNuMuBarEvents.begin();
       itBar!=fvFDMCCCNuMuBarEvents.end();
       ++itBar){
    NuEvent barEvent = *itBar;
    Int_t truthpar = this->TruthIndex(barEvent);
    Double_t barweight = 1.0;
    barweight *= barEvent.rw;
    if (truthpar>=0){
      barweight *= fvNDFitPars[truthpar];
    }
    Double_t thisSn2 = 0.0;
    Double_t thisDm2 = 0.0;
    if (1==barEvent.iaction){
      if (-14==barEvent.inu){
        thisSn2 = Sn2Bar;
        thisDm2 = Dm2Bar;
      }
      if (14==barEvent.inu){
        thisSn2 = Sn2;
        thisDm2 = Dm2;
      }
    }
    Double_t oscweightbar = 1.0-thisSn2*
      TMath::Sin(1.27*thisDm2*735.0/barEvent.neuEnMC)*
      TMath::Sin(1.27*thisDm2*735.0/barEvent.neuEnMC);
    barweight *= oscweightbar;
    Double_t recoNuE = shwEShift*(barEvent.shwEn)
      + trkEShift*(barEvent.trkEn);
    fhFDMCRecoENuMuBarCC->Fill(recoNuE,barweight);
  }
  if (ffdMCPoT){fhFDMCRecoENuMuBarCC->Scale(ffdDataPoT/ffdMCPoT);}

  //Fill the numu histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itNu
         = fvFDMCCCNuMuEvents.begin();
       itNu!=fvFDMCCCNuMuEvents.end();
       ++itNu){
    NuEvent nuEvent = *itNu;
    Int_t truthpar = this->TruthIndex(nuEvent);
    Double_t nuweight = 1.0;
    nuweight *= nuEvent.rw;
    if (truthpar>=0){
      nuweight *= fvNDFitPars[truthpar];
    }
    Double_t thisSn2 = 0.0;
    Double_t thisDm2 = 0.0;
    if (1==nuEvent.iaction){
      if (-14==nuEvent.inu){
        thisSn2 = Sn2Bar;
        thisDm2 = Dm2Bar;
      }
      if (14==nuEvent.inu){
        thisSn2 = Sn2;
        thisDm2 = Dm2;
      }
    }
    Double_t oscweight = 1.0-thisSn2*
      TMath::Sin(1.27*thisDm2*735.0/nuEvent.neuEnMC)*
      TMath::Sin(1.27*thisDm2*735.0/nuEvent.neuEnMC);
    nuweight *= oscweight;
    Double_t recoNuE = shwEShift*(nuEvent.shwEn)
      + trkEShift*(nuEvent.trkEn);
    fhFDMCRecoENuMuCC->Fill(recoNuE,nuweight);
  }
  if (ffdMCPoT){fhFDMCRecoENuMuCC->Scale(ffdDataPoT/ffdMCPoT);}
}

//____________________________________________________________________72
void NuXFitAnalysis::FillTransitionFDPrediction
(const vector<double>& par) const
{
  //Reset recoE histograms
  fhFDMCRecoENuMuCC->Reset();
  fhFDMCRecoENuMuBarCC->Reset();

  Double_t shwEShift=1.0;
  Double_t trkEShift = 1.0;

  Double_t Dm2 = fdm2;
  Double_t Sn2 = fsn2;
  Double_t Dm2Bar = fdm2;
  Double_t Sn2Bar = fsn2;
  Double_t transitionProb = par[0];

  //Fill the nubar histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itBar
         = fvFDMCCCNuMuBarEvents.begin();
       itBar!=fvFDMCCCNuMuBarEvents.end();
       ++itBar){
    NuEvent barEvent = *itBar;
    Int_t truthpar = this->TruthIndex(barEvent);
    Double_t barweight = 1.0;
    barweight *= barEvent.rw;
    if (truthpar>=0){
      barweight *= fvNDFitPars[truthpar];
    }
    Double_t thisSn2 = 0.0;
    Double_t thisDm2 = 0.0;
    if (1==barEvent.iaction){
      if (-14==barEvent.inu){
        thisSn2 = Sn2Bar;
        thisDm2 = Dm2Bar;
      }
      if (14==barEvent.inu){
        thisSn2 = Sn2;
        thisDm2 = Dm2;
      }
    }
    Double_t oscweightbar = 1.0-thisSn2*
      TMath::Sin(1.27*thisDm2*735.0/barEvent.neuEnMC)*
      TMath::Sin(1.27*thisDm2*735.0/barEvent.neuEnMC);
    Double_t transitionWeight = barweight;
    barweight *= oscweightbar;
    Double_t recoNuE = shwEShift*(barEvent.shwEn)
      + trkEShift*(barEvent.trkEn);
    Double_t trueNuE = barEvent.neuEnMC;
    fhFDMCRecoENuMuBarCC->Fill(recoNuE,barweight);
    if (1==barEvent.iaction){
      if (-14==barEvent.inu){
        transitionWeight *= 1.0-oscweightbar;
        transitionWeight *= transitionProb;
        Double_t nuBarSelEff = this->NuBarCCSelectionEfficiency(trueNuE);
        Double_t nuBarXSec = this->NuBarCCCrossSection(trueNuE);
        if (nuBarSelEff && nuBarXSec){
          transitionWeight /= nuBarSelEff;
          transitionWeight /= nuBarXSec;
        }
        transitionWeight *= this->NuMuCCSelectionEfficiency(trueNuE);
        transitionWeight *= this->NuMuCCCrossSection(trueNuE);
        fhFDMCRecoENuMuCC->Fill(recoNuE,transitionWeight);
      }
    }
  }
  if (ffdMCPoT){fhFDMCRecoENuMuBarCC->Scale(ffdDataPoT/ffdMCPoT);}

  //Fill the numu histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itNu
         = fvFDMCCCNuMuEvents.begin();
       itNu!=fvFDMCCCNuMuEvents.end();
       ++itNu){
    NuEvent nuEvent = *itNu;
    Int_t truthpar = this->TruthIndex(nuEvent);
    Double_t nuweight = 1.0;
    nuweight *= nuEvent.rw;
    if (truthpar>=0){
      nuweight *= fvNDFitPars[truthpar];
    }
    Double_t thisSn2 = 0.0;
    Double_t thisDm2 = 0.0;
    if (1==nuEvent.iaction){
      if (-14==nuEvent.inu){
        thisSn2 = Sn2Bar;
        thisDm2 = Dm2Bar;
      }
      if (14==nuEvent.inu){
        thisSn2 = Sn2;
        thisDm2 = Dm2;
      }
    }
    Double_t oscweight = 1.0-thisSn2*
      TMath::Sin(1.27*thisDm2*735.0/nuEvent.neuEnMC)*
      TMath::Sin(1.27*thisDm2*735.0/nuEvent.neuEnMC);
    Double_t transitionWeight = nuweight;
    nuweight *= oscweight;
    Double_t recoNuE = shwEShift*(nuEvent.shwEn)
      + trkEShift*(nuEvent.trkEn);
    Double_t trueNuE = nuEvent.neuEnMC;
    fhFDMCRecoENuMuCC->Fill(recoNuE,nuweight);
    if (1==nuEvent.iaction){
      if (14==nuEvent.inu){
        transitionWeight *= 1.0-oscweight;
        transitionWeight *= transitionProb;
        transitionWeight *= this->NuBarCCSelectionEfficiency(trueNuE);
        transitionWeight *= this->NuBarCCCrossSection(trueNuE);
        Double_t nuMuSelEff = this->NuMuCCSelectionEfficiency(trueNuE);
        Double_t nuMuXSec = this->NuMuCCCrossSection(trueNuE);
        if (nuMuSelEff && nuMuXSec){
          transitionWeight /= nuMuSelEff;
          transitionWeight /= nuMuXSec;
        }
        fhFDMCRecoENuMuBarCC->Fill(recoNuE,transitionWeight);
      }
    }
  }
  if (ffdMCPoT){fhFDMCRecoENuMuCC->Scale(ffdDataPoT/ffdMCPoT);}
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis
::NuBarCCCrossSection(const Double_t trueNuE) const
{
  //  return fNuBarCCCrossSectionGraph->Eval(trueNuE,0,"");
  return fhNuBarCCCrossSections->
    GetBinContent(fhNuBarCCCrossSections->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis
::NuBarCCSelectionEfficiency(const Double_t trueNuE) const
{
  return fhNuBarCCSelectionEfficiency->
    GetBinContent(fhNuBarCCSelectionEfficiency->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis
::NuMuCCCrossSection(const Double_t trueNuE) const
{
  //  return fNuMuCCCrossSectionGraph->Eval(trueNuE,0,"");
  return fhNuMuCCCrossSections->
    GetBinContent(fhNuMuCCCrossSections->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuXFitAnalysis
::NuMuCCSelectionEfficiency(const Double_t trueNuE) const
{
  return fhNuMuCCSelectionEfficiency->
    GetBinContent(fhNuMuCCSelectionEfficiency->FindBin(trueNuE));
}

//____________________________________________________________________72
const Int_t NuXFitAnalysis::TruthIndex(const NuEvent& info) const
{
  Int_t numnubins = fvTrueEBinsNuMuCC.size();
  Int_t numbarbins = fvTrueEBinsNuMuBarCC.size();

  Int_t type = info.inu;
  const vector<Double_t>* truthbins;
  Int_t offset = 0;
  Int_t numbins = 0;
  if (14 == type){
    truthbins = &fvTrueEBinsNuMuCC;
    numbins = numnubins;
    offset = 0;
  }
  else if (-14 == type){
    truthbins = &fvTrueEBinsNuMuBarCC;
    numbins = numbarbins;
    offset = numnubins;
  }
  else{
    return numnubins + numbarbins;
  }
  for (Int_t i=1; i<numbins; ++i){
    if ((*truthbins)[i]>info.neuEnMC){
      return i - 1 + offset;
    }
  }
  if ((14 == type) || (-14 == type)){
    return (numbins + offset);
  }
  return numnubins + numbarbins;
}

//____________________________________________________________________72
void NuXFitAnalysis::FDData(const string fdDataFileName)
{
  if (ftFDDataTree){delete ftFDDataTree; ftFDDataTree = 0;}
  ftFDDataTree = new NuTreeWrapper(fdDataFileName);

  //Get the total PoT
  ffdDataPoT = ftFDDataTree->GetPoT();
  
  MSG("NuXFitAnalysis.cxx",Msg::kInfo)
    << "fdDataPoT: " << ffdDataPoT
    << endl;
}

//____________________________________________________________________72
void NuXFitAnalysis::FDNuMuBarCCData(const TH1D& fdData)
{
  if (fhFDNuMuBarCCDataFromFile){
    delete fhFDNuMuBarCCDataFromFile; fhFDNuMuBarCCDataFromFile = 0;
  }
  fhFDNuMuBarCCDataFromFile = new TH1D(fdData);
  fhFDNuMuBarCCDataFromFile->SetDirectory(0);
}

//____________________________________________________________________72
void NuXFitAnalysis::FDNuMuCCData(const TH1D& fdData)
{
  if (fhFDNuMuCCDataFromFile){
    delete fhFDNuMuCCDataFromFile; fhFDNuMuCCDataFromFile = 0;
  }
  fhFDNuMuCCDataFromFile = new TH1D(fdData);
  fhFDNuMuCCDataFromFile->SetDirectory(0);
}

//____________________________________________________________________72
void NuXFitAnalysis::FDMC(const string fdMCFileName)
{
  if (ftFDMCTree){delete ftFDMCTree; ftFDMCTree = 0;}
  ftFDMCTree = new NuTreeWrapper(fdMCFileName);

  //Get the total PoT
  ffdMCPoT = ftFDMCTree->GetPoT();
  
  MSG("NuXFitAnalysis.cxx",Msg::kInfo)
    << "fdMCPoT: " << ffdMCPoT
    << endl;
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CreateHistograms()
{
  Bool_t goodConfig = true;
  if (!this->CreateFDNuMuCCDataHisto()) goodConfig = false;
  if (!this->CreateFDNuMuBarCCDataHisto()) goodConfig = false;
  if (!this->CreateFDNuMuCCMCHisto()) goodConfig = false;
  if (!this->CreateFDNuMuBarCCMCHisto()) goodConfig = false;
  return goodConfig;
}

//____________________________________________________________________72
void NuXFitAnalysis::ConfigureForExternalFit
(const NuXFit::NuXAnalysis_t analSetting)
{
  fanalysisSetting = analSetting;
  this->CreateHistograms();
  this->FillFDDataHistograms();
  this->CacheFDMCEvents();

  this->SetupPars();
}

//____________________________________________________________________72
void NuXFitAnalysis::CPTAnalysis()
{
  fanalysisSetting = NuXFit::kCPT;

  if (!this->CreateHistograms()){return;}
  if (!this->FillFDDataHistograms()){return;}
  if (!this->CacheFDMCEvents()){return;}

  if (!fdoGridSearch){
    this->SetupPars();
    
    fFitter->FixParameter(0);
    fFitter->FixParameter(1);
    
    fFitter->Minimize();
    fFitter->PrintResults(0,0);
    
    fBestDm2 = fFitter->GetParameter(0);
    fBestSn2 = fFitter->GetParameter(1);
    fBestDm2Bar = fFitter->GetParameter(2);
    fBestSn2Bar = fFitter->GetParameter(3);
  }

  if (fdoGridSearch){
    Double_t bestChi2 = 1.0e10;
    Double_t bestsn2bar = 0.0;
    Double_t bestdm2bar = 0.0;
    
    Double_t sn2BarLow = 0.8;
    Double_t sn2BarHigh = 1.0;
    Double_t sn2BarGranularity = 0.001;
    Double_t dm2BarLow = 1.0e-3;
    Double_t dm2BarHigh = 4.0e-3;
    Double_t dm2BarGranularity = 0.005e-3; 
    
    //Create the Chi2 histogram
    Int_t numSn2BarBins = (Int_t) round((sn2BarHigh-sn2BarLow)/sn2BarGranularity) + 1;
    Int_t numDm2BarBins = (Int_t) round((dm2BarHigh-dm2BarLow)/dm2BarGranularity) + 1;
    TH2D hChi2Bar("hChi2Bar",
                  "",
                  numSn2BarBins,
                  sn2BarLow-sn2BarGranularity/2.0,
                  sn2BarHigh+sn2BarGranularity/2.0,
                  numDm2BarBins,
                  dm2BarLow-dm2BarGranularity/2.0,
                  dm2BarHigh+dm2BarGranularity/2.0); 
    
    for (Double_t sn2bar = sn2BarLow;
         sn2bar <= sn2BarHigh+sn2BarGranularity/2.0;
         sn2bar += sn2BarGranularity){
      MSG("NuMatrixFitter.cxx",Msg::kInfo) << "sn2bar = " << sn2bar << endl;
      for (Double_t dm2bar = dm2BarLow;
           dm2bar <= dm2BarHigh+dm2BarGranularity/2.0;
           dm2bar += dm2BarGranularity){
        
        vector<Double_t> vParameters;
        vParameters.push_back(2.38);
        vParameters.push_back(1.0);
        vParameters.push_back(dm2bar);
        vParameters.push_back(sn2bar);
        this->FillCPTFDPrediction(vParameters);
        
        Double_t chi2 = this->CalculateLikelihood(fhFDMCRecoENuMuBarCC,
                                                  fhFDDataRecoENuMuBarCC);
        
        //       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

    fBestDm2 = 2.38;
    fBestSn2 = 1.0;
    fBestDm2Bar = bestdm2bar;
    fBestSn2Bar = bestsn2bar;
  }

  vector<Double_t> vBestFit;
  vBestFit.push_back(fBestDm2);
  vBestFit.push_back(fBestSn2);
  vBestFit.push_back(fBestDm2Bar);
  vBestFit.push_back(fBestSn2Bar);
  this->FillCPTFDPrediction(vBestFit);
  if (fwriteOutput){
    fOutputFile->cd();
    fhFDMCRecoENuMuCC->Write("BestFitFDMCRecoENuMuCC");
    fhFDMCRecoENuMuBarCC->Write("BestFitFDMCRecoENuMuBarCC");
  }



  
//   ROOT::Minuit2::MnUserParameters mnPars;
//   mnPars.Add("Dm2",0.003,0.2,0.0e-3,10.0e-3);
//   mnPars.Add("Sn2",1.0,0.2,0.0,2.0);
//   mnPars.Add("Dm2Bar",0.003,0.2,4.0e-3,8.0e-3);
//   mnPars.Add("Sn2Bar",1.0,0.2,0.4,1.6);
//   ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
//   cout << "Beginning second fit" << endl;
//   ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();
//   cout << "Function minimum is " << funcMin.UserParameters().Parameters()[2].Value() << ", " << funcMin.UserParameters().Parameters()[3].Value() << endl;
  
//   ROOT::Minuit2::FCNBase& fcnBase = *this;
//   ROOT::Minuit2::MnContours contours(fcnBase,funcMin);
//   this->SetErrorDef(4);
//   ContoursError contErr = contours.Contour(2,3,20);
//   //  cout << contErr << endl;
//   cout << "Contour:" << endl;

//   std::vector<std::pair<Double_t,Double_t> > vContPoints = contErr();

//   Double_t* x = new Double_t[vContPoints.size()];
//   for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
//     x[pointCount] = vContPoints[pointCount].first;
//   }
//   Double_t* y = new Double_t[vContPoints.size()];
//   for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
//     y[pointCount] = vContPoints[pointCount].second;
//   }
//   TGraph* gCont = new TGraph(vContPoints.size(),y,x);
//   TCanvas* cCont;
//   gCont->Draw("AL");
  //  cout << contours << endl;

  /*
  TGraph2DErrors* gConf = new TGraph2DErrors(100);
  fFitter->FixParameter(0);
  fFitter->FixParameter(1);
  fFitter->GetFitter()->GetConfidenceIntervals(gConf,0.66);
  if (fwriteOutput){
    fOutputFile->cd();
    gConf->Write("OneSigmaContour");
  }
  */
}

//____________________________________________________________________72
void NuXFitAnalysis::TransitionAnalysis(const Double_t dm2,
                                        const Double_t sn2)
{
  fanalysisSetting = NuXFit::kTransition;
  fdm2 = dm2;
  fsn2 = sn2;

  if (!this->CacheFDMCEvents()){return;}
  if (!this->FillFDDataHistograms()){return;}

  if (!fdoGridSearch){
    this->SetupPars();
    
    fFitter->Minimize();
    fFitter->PrintResults(0,0);
    
    fBestTransitionProb = fFitter->GetParameter(0);
  }

  if (fdoGridSearch){
    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((transitionProbHigh-transitionProbLow)/
                    transitionProbGranularity) + 1;
    TH1D hChi2Trans("hChi2Trans",
                    "",
                    numTransitionProbBins,
                    transitionProbLow-transitionProbGranularity/2.0,
                    transitionProbHigh+transitionProbGranularity/2.0);
    for (Double_t transitionProb = transitionProbLow;
         transitionProb <= transitionProbHigh+transitionProbGranularity/2.0;
         transitionProb += transitionProbGranularity){
      MSG("NuMatrixFitter.cxx",Msg::kInfo)
        << "Fitting transitionProb = " << transitionProb << endl;

      vector<Double_t> vParams;
      vParams.push_back(transitionProb);
      this->FillTransitionFDPrediction(vParams);
      
      Double_t chi2 = this->CalculateLikelihood(fhFDMCRecoENuMuBarCC,
                                                fhFDDataRecoENuMuBarCC);
      
      Int_t xBin = hChi2Trans.GetXaxis()->FindBin(transitionProb);
      hChi2Trans.SetBinContent(xBin,chi2);
      if (chi2 < bestChi2){
        bestChi2 = chi2;
        bestTransitionProb = transitionProb;
      }
    } //transitionProb
    fBestTransitionProb = bestTransitionProb;
  }
  
  vector<Double_t> vBestFit;
  vBestFit.push_back(fBestTransitionProb);
  this->FillTransitionFDPrediction(vBestFit);
  if (fwriteOutput){
    fOutputFile->cd();
    fhFDMCRecoENuMuCC->Write("BestFitFDMCRecoENuMuCC");
    fhFDMCRecoENuMuBarCC->Write("BestFitFDMCRecoENuMuBarCC");
  }
}

//____________________________________________________________________72
void NuXFitAnalysis::SetupPars()
{
  if (NuXFit::kCPT == fanalysisSetting){
    this->SetupCPTPars();
    return;
  }
  else if (NuXFit::kTransition == fanalysisSetting){
    this->SetupTransitionPars();
    return;
  }
  else{
    MSG("NuXFitAnalysis.cxx",Msg::kFatal)
      << "Incorrect analysis setting" << endl;
    return;
  }
}

//____________________________________________________________________72
void NuXFitAnalysis::SetupCPTPars()
{
  fFitter->SetParameter(0,
                        "Dm2",
                        0.00238,0.2,1.0,0.0);
  fFitter->SetParameter(1,
                        "Sn2",
                        1.0,0.2,1.0,0.0);
  fFitter->SetParameter(2,
                        "Dm2Bar",
                        0.00238,0.2,1.0,0.0);
  fFitter->SetParameter(3,
                        "Sn2Bar",
                        1.0,0.2,1.0,0.0);
}

//____________________________________________________________________72
void NuXFitAnalysis::SetupTransitionPars()
{
  //Get cross-sections
  TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");

  if (fNuMuCCCrossSectionGraph){
    delete fNuMuCCCrossSectionGraph;
    fNuMuCCCrossSectionGraph = 0;
  }
  if (fhNuMuCCCrossSections){
    delete fhNuMuCCCrossSections;
    fhNuMuCCCrossSections = 0;
  }
  TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
  fhNuMuCCCrossSections = new TH1F(*XSec_CC);
  fhNuMuCCCrossSections->SetDirectory(0);
  Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
  Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
  for(int i=0;i<XSec_CC->GetNbinsX();i++) {
    x[i] = XSec_CC->GetBinCenter(i+1);
    y[i] = XSec_CC->GetBinContent(i+1);
  }
  fNuMuCCCrossSectionGraph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}

  if (fNuBarCCCrossSectionGraph){
    delete fNuBarCCCrossSectionGraph;
    fNuBarCCCrossSectionGraph = 0;
  }
  if (fhNuBarCCCrossSections){
    delete fhNuBarCCCrossSections;
    fhNuBarCCCrossSections = 0;
  }
  XSec_CC = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
  fhNuBarCCCrossSections = new TH1F(*XSec_CC);
  fhNuBarCCCrossSections->SetDirectory(0);
  x = new Float_t[XSec_CC->GetNbinsX()];
  y = new Float_t[XSec_CC->GetNbinsX()];
  for(int i=0;i<XSec_CC->GetNbinsX();i++) {
    x[i] = XSec_CC->GetBinCenter(i+1);
    y[i] = XSec_CC->GetBinContent(i+1);
  }
  fNuBarCCCrossSectionGraph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}
  
  xsecfile->Close();
  if (xsecfile){delete xsecfile; xsecfile = 0;}
  
  //Get selection efficiencies
  TFile* selEffFile = new TFile(fSelEffFileName.c_str(),"READ");

  if (fhNuBarCCSelectionEfficiency){
    delete fhNuBarCCSelectionEfficiency;
    fhNuBarCCSelectionEfficiency = 0;
  }
  fhNuBarCCSelectionEfficiency = (TH1D*) selEffFile->Get("EfficiencyPQ_FD");
  fhNuBarCCSelectionEfficiency->SetDirectory(0);

  if (fhNuMuCCSelectionEfficiency){
    delete fhNuMuCCSelectionEfficiency;
    fhNuMuCCSelectionEfficiency = 0;
  }
  fhNuMuCCSelectionEfficiency = (TH1D*) selEffFile->Get("Efficiency_FD");
  fhNuMuCCSelectionEfficiency->SetDirectory(0);

  selEffFile->Close();
  if (selEffFile){delete selEffFile; selEffFile = 0;}

  //Setup fit parameters
  fFitter->SetParameter(0,
                        "transitionProb",
                        0.0,0.2,1.0,0.0);
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::FillFDDataHistograms()
{
  //Preferentially take data from histograms if they've been supplied
  if (fhFDNuMuCCDataFromFile && fhFDNuMuBarCCDataFromFile){
    if (fhFDDataRecoENuMuCC){
      delete fhFDDataRecoENuMuCC; fhFDDataRecoENuMuCC = 0;
    }
    if (fhFDDataRecoENuMuBarCC){
      delete fhFDDataRecoENuMuBarCC; fhFDDataRecoENuMuBarCC = 0;
    }
    fhFDDataRecoENuMuCC = new TH1D(*fhFDNuMuCCDataFromFile);
    fhFDDataRecoENuMuBarCC = new TH1D(*fhFDNuMuBarCCDataFromFile);

    if (fwriteOutput){
      fOutputFile->cd();
      fhFDDataRecoENuMuCC->Write("InputFDNuMuCCData");
      fhFDDataRecoENuMuBarCC->Write("InputFDNuMuBarCCData");
    }
    return true;
  }

  if (!fhFDDataRecoENuMuCC){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD numu data histogram not created"
      << endl;
    return false;
  }
  if (!fhFDDataRecoENuMuBarCC){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD numubar data histogram not created"
      << endl;
    return false;
  }
  if (!ftFDDataTree){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD data not supplied"
      << endl;
    return false;
  }

  fhFDDataRecoENuMuCC->Reset();
  fhFDDataRecoENuMuBarCC->Reset();
  
  Double_t Sn2 = 1.0;
  Double_t Sn2Bar = 1.0;
  Double_t Dm2 = 0.003;
  Double_t Dm2Bar = 0.003;
  /*
  Double_t Sn2 = 0.0;
  Double_t Sn2Bar = 0.0;
  Double_t Dm2 = 0.0;
  Double_t Dm2Bar = 0.0;
  */
  MSG("NuXFitAnalysis.cxx",Msg::kInfo)
    << "Looping over FD data:"
    << endl;
  Int_t dataEntries = ftFDDataTree->GetEntries();
  for (Int_t counter = 0;
       counter < dataEntries;
       ++counter){
    if (!(counter%10000)){
      MSG("NuXFitAnalysis.cxx",Msg::kInfo)
        << counter << " of " << dataEntries
        << endl;
    }
    NuEvent currentEvent = ftFDDataTree->GetInfoObject(counter);
    NuSelectionResult_t recoType = this->SelectedAs(currentEvent);
    Double_t recoEnergy = currentEvent.energy;
    Double_t nuweight = 1.0;
    nuweight *= currentEvent.rw;
    
    //Fake oscillations
    Double_t thisSn2 = 0.0;
    Double_t thisDm2 = 0.0;
    if (1==currentEvent.iaction){
      if (-14==currentEvent.inu){
        thisSn2 = Sn2Bar;
        thisDm2 = Dm2Bar;
      }
      if (14==currentEvent.inu){
        thisSn2 = Sn2;
        thisDm2 = Dm2;
      }
    }
    Double_t oscweight = 1.0-thisSn2*
      TMath::Sin(1.27*thisDm2*735.0/currentEvent.neuEnMC)*
      TMath::Sin(1.27*thisDm2*735.0/currentEvent.neuEnMC);
    nuweight *= oscweight;
    
    if (recoType==NuXSecFit::kCCNuMu){
      fhFDDataRecoENuMuCC->Fill(recoEnergy,nuweight);
    }
    if (recoType==NuXSecFit::kCCNuMuBar){
      fhFDDataRecoENuMuBarCC->Fill(recoEnergy,nuweight);
    }
  }

  if (fwriteOutput){
    fOutputFile->cd();
    fhFDDataRecoENuMuCC->Write("InputFDNuMuCCData");
    fhFDDataRecoENuMuBarCC->Write("InputFDNuMuBarCCData");
  }

  return true;
}

//____________________________________________________________________72
const Bool_t NuXFitAnalysis::CacheFDMCEvents()
{
  if (!ftFDMCTree){
    MSG("NuXFitAnalysis.cxx",Msg::kWarning)
      << "FD MC not supplied"
      << endl;
    return false;
  }

  TH1D hFDRawMCRecoECCNuMuBar(*fhFDDataRecoENuMuBarCC);
  TH1D hFDRawMCRecoECCNuMu(*fhFDDataRecoENuMuCC);
  TH1D hFDNoOscPredRecoECCNuMuBar(*fhFDDataRecoENuMuBarCC);
  TH1D hFDNoOscPredRecoECCNuMu(*fhFDDataRecoENuMuCC);

  fvFDMCCCNuMuEvents.clear();
  fvFDMCCCNuMuBarEvents.clear();

  MSG("NuXFitAnalysis.cxx",Msg::kInfo)
    << "Looping over FD MC:"
    << endl;
  Int_t mcEntries = ftFDMCTree->GetEntries();
  for (Int_t counter = 0;
       counter < mcEntries;
       ++counter){
    if (!(counter%10000)){
      MSG("NuXFitAnalysis.cxx",Msg::kInfo)
        << counter << " of " << mcEntries
        << endl;
    }
    NuEvent currentEvent = ftFDMCTree->GetInfoObject(counter);
    fnuSystematic->Shift(currentEvent);
    NuSelectionResult_t recoType = this->SelectedAs(currentEvent);
    Double_t eventWeight = currentEvent.rw;

    Int_t truthpar = this->TruthIndex(currentEvent);
    Double_t eventWeightPred = eventWeight;
    if (truthpar>=0){
      eventWeightPred *= fvNDFitPars[truthpar];
    }

    if (NuXSecFit::kCCNuMu==recoType){
      fvFDMCCCNuMuEvents.push_back(currentEvent);
      hFDRawMCRecoECCNuMu.Fill(currentEvent.energy,eventWeight);
      hFDNoOscPredRecoECCNuMu.Fill(currentEvent.energy,eventWeightPred);
    }
    if (NuXSecFit::kCCNuMuBar==recoType){
      fvFDMCCCNuMuBarEvents.push_back(currentEvent);
      hFDRawMCRecoECCNuMuBar.Fill(currentEvent.energy,eventWeight);
      hFDNoOscPredRecoECCNuMuBar.Fill(currentEvent.energy,eventWeightPred);
    }
  }
  
  if (ffdMCPoT){
    hFDNoOscPredRecoECCNuMu.Scale(ffdDataPoT/ffdMCPoT);
    hFDNoOscPredRecoECCNuMuBar.Scale(ffdDataPoT/ffdMCPoT);
    hFDRawMCRecoECCNuMu.Scale(ffdDataPoT/ffdMCPoT);
    hFDRawMCRecoECCNuMuBar.Scale(ffdDataPoT/ffdMCPoT);
  }
  
  if (fwriteOutput){
    fOutputFile->cd();
    hFDNoOscPredRecoECCNuMu.Write("NoOscPredRecoECCNuMu");
    hFDNoOscPredRecoECCNuMuBar.Write("NoOscPredRecoECCNuMuBar");
    hFDRawMCRecoECCNuMu.Write("RawMCRecoECCNuMu");
    hFDRawMCRecoECCNuMuBar.Write("RawMCRecoECCNuMuBar");
  }

  return true;
}

//____________________________________________________________________72
const NuSelectionResult_t
NuXFitAnalysis::SelectedAs(NuEvent& nuEvent) const
{
  NuCuts nuCuts;
  nuEvent.anaVersion = this->SelectionScheme();
  
  //cut on LI
  if (nuCuts.IsLI(nuEvent)) return NuXSecFit::kUnknown;
  //cut on the data quality
  if (!nuCuts.IsGoodDataQuality(nuEvent)) return NuXSecFit::kUnknown;
  //cut on the spill time
  if (!nuCuts.IsGoodTimeToNearestSpill(nuEvent)) return NuXSecFit::kUnknown;
  //cut on the beam
  if (!nuCuts.IsGoodBeam(nuEvent)) return NuXSecFit::kUnknown;
  
  //cut on whether event vtx is in fid vol (does nothing for CC ana)
  //    if (!nuCuts.IsInFidVol(evt,nu)) continue;
  //ensure good number of tracks in event (note preselection above)
  if (!nuCuts.IsGoodNumberOfTracks(nuEvent)) return NuXSecFit::kUnknown;
  //check trk is in fiducial volume (note preselection above)
  if (!nuCuts.IsInFidVolTrk(nuEvent)) return NuXSecFit::kUnknown;
  
  //require a good trk fit
  if (!nuCuts.IsGoodTrackFitPass(nuEvent)) return NuXSecFit::kUnknown;
  //require a forward going neutrino about beam direction
  if (!nuCuts.IsGoodDirCos(nuEvent)) return NuXSecFit::kUnknown;
  
  //cut on the PID
  if (!nuCuts.IsGoodPID(nuEvent)) {
    return NuXSecFit::kUnknown;
  }
  //cut on the fractional track momentum and sign error      
  if (!nuCuts.IsGoodSigmaQP_QP(nuEvent)) {
    return NuXSecFit::kUnknown;
  }
  //cut on the track fit probability      
  if (!nuCuts.IsGoodFitProb(nuEvent)) {
    return NuXSecFit::kUnknown;
  }
  if (!nuCuts.IsGoodMajorityCurvature(nuEvent)){
    return NuXSecFit::kUnknown;
  }
  
  if (-1==nuEvent.charge){return NuXSecFit::kCCNuMu;}
  if (1==nuEvent.charge){return NuXSecFit::kCCNuMuBar;}
  return NuXSecFit::kUnknown;
}

---