NuCrossSectionFitter.cxx

Go to the documentation of this file.

// $Id: NuCrossSectionFitter.cxx,v 1.10 2008/05/02 19:36:16 evans Exp $
//
// Combined numu & numubar extrapolation by cross-section fitting
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
//
// $Log: NuCrossSectionFitter.cxx,v $
// Revision 1.10  2008/05/02 19:36:16  evans
// Correcting a comment.
//
// Revision 1.9  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.8  2008/01/10 17:58:33  evans
// Adding IsGoodMajorityCurvature to the list of selection cuts.
//
// Revision 1.7  2007/12/21 17:41:18  evans
// Debugging and creating a FD output file for the XFit CPT analysis.
//
// Revision 1.6  2007/12/20 19:12:46  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.5  2007/12/07 14:53:15  evans
// Adding protection against uninitialised variables found by Robert
// Hatcher when compiling with gcc+optimization.
//
// Revision 1.4  2007/12/04 12:03:07  evans
// Removing double construnction of NuSystematic object.
//
// Revision 1.3  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.2  2007/11/15 13:58:15  evans
// Removing three unused parameters.
//
// Revision 1.1  2007/11/15 13:45:25  evans
// The class to do the NDFit part of my cross-section fitting extrapolation.
//
//

#include "Minuit2/MinuitParameter.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TMinuit.h"

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

ClassImp(NuCrossSectionFitter)

CVSID("$Id: NuCrossSectionFitter.cxx,v 1.10 2008/05/02 19:36:16 evans Exp $");

using namespace NuXSecFit;

using std::string;
using std::vector;
using ROOT::Minuit2::MinuitParameter;
using ROOT::Minuit2::MnUserParameterState;

//____________________________________________________________________72
NuCrossSectionFitter::NuCrossSectionFitter()
  : fndDataPoT(0.0),
    fndMCPoT(0.0)
{
//   fanaVersion = NuCuts::kCC0250Std;
  fanaVersion = NuCuts::kJJE1;

  fnuSystematic = new NuSystematic();

  fhNDMCRecoENuMuCC = 0;
  fhNDMCRecoENuMuBarCC = 0;
  fhNDDataRecoENuMuCC = 0;
  fhNDDataRecoENuMuBarCC = 0;
  fhNDRawMCTrueECCNuMu = 0;
  fhNDRawMCTrueECCNuMuBar = 0;
  fhNDRawMCRecoECCNuMu = 0;
  fhNDRawMCRecoECCNuMuBar = 0;

  fhNDNuMuCCDataFromFile = 0;
  fhNDNuMuBarCCDataFromFile = 0;

  fOutputFile = 0;

  ftNDDataTree = 0;
  ftNDMCTree = 0;

  fFitter = new TFitterMinuit();
  fFitter->CreateMinimizer();
  fFitter->SetMinuitFCN(this);
  
  //ND fake parameters (set to 1.0)
  fvNDFakeTrueENuMuCCBins.push_back(30.0);
  fvNDFakeTrueENuMuBarCCBins.push_back(30.0);
  for (Int_t counter=0; counter<10; ++counter){
    fvfakepars.push_back(1.0);
  }
}

//____________________________________________________________________72
NuCrossSectionFitter::NuCrossSectionFitter(const NuXMLConfig& xmlConfig)
  : fndDataPoT(0.0),
    fndMCPoT(0.0)
{
//   fanaVersion = NuCuts::kCC0250Std;
  fanaVersion = NuCuts::kJJE1;

  fnuSystematic = new NuSystematic(xmlConfig);

  fhNDMCRecoENuMuCC = 0;
  fhNDMCRecoENuMuBarCC = 0;
  fhNDDataRecoENuMuCC = 0;
  fhNDDataRecoENuMuBarCC = 0;
  fhNDRawMCTrueECCNuMu = 0;
  fhNDRawMCTrueECCNuMuBar = 0;
  fhNDRawMCRecoECCNuMu = 0;
  fhNDRawMCRecoECCNuMuBar = 0;

  fhNDNuMuCCDataFromFile = 0;
  fhNDNuMuBarCCDataFromFile = 0;

  fOutputFile = 0;

  ftNDDataTree = 0;
  ftNDMCTree = 0;

  fFitter = new TFitterMinuit();
  fFitter->CreateMinimizer();
  fFitter->SetMinuitFCN(this);
  
  //ND fake parameters (set to 1.0)
  fvNDFakeTrueENuMuCCBins.push_back(30.0);
  fvNDFakeTrueENuMuBarCCBins.push_back(30.0);
  for (Int_t counter=0; counter<10; ++counter){
    fvfakepars.push_back(1.0);
  }
}

//____________________________________________________________________72
NuCrossSectionFitter::~NuCrossSectionFitter()
{
  if (fOutputFile){
    fOutputFile->Close();
    delete fOutputFile;
    fOutputFile = 0;
  }
  if (fhNDMCRecoENuMuCC){delete fhNDMCRecoENuMuCC; fhNDMCRecoENuMuCC=0;}
  if (fhNDMCRecoENuMuBarCC){
    delete fhNDMCRecoENuMuBarCC; fhNDMCRecoENuMuBarCC = 0;
  }
  if (fhNDDataRecoENuMuCC){
    delete fhNDDataRecoENuMuCC; fhNDDataRecoENuMuCC = 0;
  }
  if (fhNDDataRecoENuMuBarCC){
    delete fhNDDataRecoENuMuBarCC; fhNDDataRecoENuMuBarCC = 0;
  }
  if (fhNDRawMCTrueECCNuMu){
    delete fhNDRawMCTrueECCNuMu; fhNDRawMCTrueECCNuMu = 0;
  }
  if (fhNDRawMCTrueECCNuMuBar){
    delete fhNDRawMCTrueECCNuMuBar; fhNDRawMCTrueECCNuMuBar = 0;
  }
  if (fhNDRawMCRecoECCNuMu){
    delete fhNDRawMCRecoECCNuMu; fhNDRawMCRecoECCNuMu = 0;
  }
  if (fhNDRawMCRecoECCNuMuBar){
    delete fhNDRawMCRecoECCNuMuBar; fhNDRawMCRecoECCNuMuBar = 0;
  }

  if (fhNDNuMuCCDataFromFile){
    delete fhNDNuMuCCDataFromFile; fhNDNuMuCCDataFromFile = 0;
  }
  if (fhNDNuMuBarCCDataFromFile){
    delete fhNDNuMuBarCCDataFromFile; fhNDNuMuBarCCDataFromFile = 0;
  }

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

  //Note that the object fFitter newed in the constructor is never
  //deleted. This is an intentional memory leak to hack my way around
  //typical ROOT stupidity. I have to pass fFitter a this pointer to
  //tell it to use this object for its minimisation. If I delete
  //fFitter, it promptly runs the destructor of this. This causes an
  //infinite chain of destructors which, unsurprisingy, causes a seg
  //fault.
}

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

//____________________________________________________________________72
double NuCrossSectionFitter::operator () (const vector<double>& pars) const
{
  static Int_t fitcounter = 0;
  if (!(fitcounter%50)){
    cout << "Minuit TF call " << fitcounter << endl;
  }
  ++fitcounter;
  this->FillTruthWeightedHistograms(pars);
  return this->CalculateLikelihood();
}

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

  Double_t like = 0;

  for (Int_t i=1; i<=fhNDMCRecoENuMuCC->GetNbinsX(); ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fhNDMCRecoENuMuCC->GetBinContent(i);
    Double_t dnu = fhNDDataRecoENuMuCC->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<=fhNDMCRecoENuMuBarCC->GetNbinsX(); ++i){
    Double_t mbar = fhNDMCRecoENuMuBarCC->GetBinContent(i);
    Double_t dbar = fhNDDataRecoENuMuBarCC->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 NuCrossSectionFitter::FillTruthWeightedHistograms
(const vector<double>& par) const
{
  //Reset recoE histograms (removes bin contents and errors)
  fhNDMCRecoENuMuCC->Reset();
  fhNDMCRecoENuMuBarCC->Reset();

  Int_t numnubins = fvNDTrueENuMuCCBins.size();
  Int_t numbarbins = fvNDTrueENuMuBarCCBins.size();
  //NDFit parameters:
  Int_t shwEpar = numnubins + numbarbins + 1;
  Double_t shwEShift = par[shwEpar];
  Int_t trkEpar = numnubins + numbarbins + 2;
  Double_t trkEShift = par[trkEpar];

  //Fill the nubar histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itBar
         = fvNDMCCCNuMuBarEvents.begin();
       itBar!=fvNDMCCCNuMuBarEvents.end();
       ++itBar){
    NuEvent barEvent = *itBar;
    Int_t truthpar = this->TruthIndex(barEvent);
    Double_t barweight = 1.0;
    barweight *= barEvent.rw;
    if (truthpar>=0){
      barweight *= par[truthpar];
    }

    Double_t recoNuE = shwEShift*(barEvent.shwEn)
      + trkEShift*(barEvent.trkEn);
    fhNDMCRecoENuMuBarCC->Fill(recoNuE,barweight);
  }
  if (fndMCPoT){fhNDMCRecoENuMuBarCC->Scale(fndDataPoT/fndMCPoT);}

  //Fill the numu histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itNu
         = fvNDMCCCNuMuEvents.begin();
       itNu!=fvNDMCCCNuMuEvents.end();
       ++itNu){
    NuEvent nuEvent = *itNu;
    Int_t truthpar = this->TruthIndex(nuEvent);
    Double_t nuweight = 1.0;
    nuweight *= nuEvent.rw;
    if (truthpar>=0){
      nuweight *= par[truthpar];
    }

    Double_t recoNuE = shwEShift*(nuEvent.shwEn)
      + trkEShift*(nuEvent.trkEn);
    fhNDMCRecoENuMuCC->Fill(recoNuE,nuweight);
  }
  if (fndMCPoT){fhNDMCRecoENuMuCC->Scale(fndDataPoT/fndMCPoT);}
}

//____________________________________________________________________72
void NuCrossSectionFitter::FillTruthWeightedHistogramsSmooth
(const vector<double>& par) const
{
  //Reset recoE histograms
  fhNDMCRecoENuMuCC->Reset();
  fhNDMCRecoENuMuBarCC->Reset();

  Int_t numnubins = fvNDTrueENuMuCCBins.size();
  Int_t numbarbins = fvNDTrueENuMuBarCCBins.size();
  //NDFit parameters:
  Int_t shwEpar = numnubins + numbarbins + 1;
  Double_t shwEShift = par[shwEpar];
  Int_t trkEpar = numnubins + numbarbins + 2;
  Double_t trkEShift = par[trkEpar];

  //Fill the nubar histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itBar
         = fvNDMCCCNuMuBarEvents.begin();
       itBar!=fvNDMCCCNuMuBarEvents.end();
       ++itBar){
    NuEvent barEvent = *itBar;
    Int_t truthpar = this->TruthIndex(barEvent);
    Double_t barweight = 1.0;
    barweight *= barEvent.rw;
    if (truthpar>=0){
      barweight *= par[truthpar];
    }

    Double_t recoNuE = shwEShift*(barEvent.shwEn)
      + trkEShift*(barEvent.trkEn);



    //Which bin am I in?
    Int_t thisBin = fhNDMCRecoENuMuBarCC->GetXaxis()->FindFixBin(recoNuE);
    //What's the bin width?
    Double_t binWidth = fhNDMCRecoENuMuBarCC->GetXaxis()->GetBinWidth(thisBin);
    //Where's the bin center?
    Double_t binCenter = fhNDMCRecoENuMuBarCC->GetXaxis()->GetBinCenter(thisBin);
    //What fraction of me goes in this bin?
    Double_t fracThis = 
      1.0 - (fabs(recoNuE-binCenter))/binWidth;
    if ((fhNDMCRecoENuMuBarCC->GetXaxis()->GetNbins()+1)==thisBin){
      fracThis = 1.0;
    }
        if (fracThis<0.0){cout << "Too many bar events!!!!!!!!!!!!!!!" << fracThis << endl;}
    fhNDMCRecoENuMuBarCC->AddBinContent(thisBin,fracThis*barweight);
    //Am I above or below bin center?
    if (recoNuE>=binCenter){
      fhNDMCRecoENuMuBarCC->AddBinContent(thisBin+1,(1.0-fracThis)*barweight);
    }
    else{
      if (thisBin>1){
        fhNDMCRecoENuMuBarCC->AddBinContent(thisBin-1,(1.0-fracThis)*barweight);
      }
      else{
        fhNDMCRecoENuMuBarCC->AddBinContent(thisBin,(1.0-fracThis)*barweight);
      }
    }



    this->FillSmooth(fhNDMCRecoENuMuBarCC,recoNuE,barweight);
  }
  if (fndMCPoT){fhNDMCRecoENuMuBarCC->Scale(fndDataPoT/fndMCPoT);}

  //Fill the numu histogram, weighted with truth parameters
  for (vector<NuEvent>::const_iterator itNu
         = fvNDMCCCNuMuEvents.begin();
       itNu!=fvNDMCCCNuMuEvents.end();
       ++itNu){
    NuEvent nuEvent = *itNu;
    Int_t truthpar = this->TruthIndex(nuEvent);
    Double_t nuweight = 1.0;
    nuweight *= nuEvent.rw;
    if (truthpar>=0){
      nuweight *= par[truthpar];
    }

    Double_t recoNuE = shwEShift*(nuEvent.shwEn)
      + trkEShift*(nuEvent.trkEn);
    fhNDMCRecoENuMuCC->Fill(recoNuE,nuweight);
  }
  if (fndMCPoT){fhNDMCRecoENuMuCC->Scale(fndDataPoT/fndMCPoT);}
}

//____________________________________________________________________72
void NuCrossSectionFitter::FillSmooth(const TH1D*/* histogram*/,
                                      Double_t /*energy*/,
                                      Double_t /*evtWeight*/) const
{
}

//____________________________________________________________________72
const Int_t NuCrossSectionFitter::TruthIndex(const NuEvent& info) const
{
  Int_t numnubins = fvNDTrueENuMuCCBins.size();
  Int_t numbarbins = fvNDTrueENuMuBarCCBins.size();

  Int_t type = info.inu;
  const vector<Double_t>* truthbins;
  Int_t offset = 0;
  Int_t numbins = 0;
  if (14 == type){
    truthbins = &fvNDTrueENuMuCCBins;
    numbins = numnubins;
    offset = 0;
  }
  else if (-14 == type){
    truthbins = &fvNDTrueENuMuBarCCBins;
    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
const Int_t NuCrossSectionFitter::FakeTruthIndex(const NuEvent& info) const
{
  Int_t numnubins = fvNDFakeTrueENuMuCCBins.size();
  Int_t numbarbins = fvNDFakeTrueENuMuBarCCBins.size();

  Int_t type = info.inu;
  const vector<Double_t>* truthbins;
  Int_t offset = 0;
  Int_t numbins = 0;
  if (14 == type){
    truthbins = &fvNDFakeTrueENuMuCCBins;
    numbins = numnubins;
    offset = 0;
  }
  else if (-14 == type){
    truthbins = &fvNDFakeTrueENuMuBarCCBins;
    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 NuCrossSectionFitter::CalculateFDWeights()
{
  const Bool_t goodConfig = this->Config();
  if (!goodConfig){
    MSG("NuCrossSectionFitter",Msg::kFatal)
      << "About to die due to bad config."
      << endl;
    return;
  }
  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "Configured successfully. Starting XFit."
    << endl;

  this->DoNDFit();

  //  this->FillFDPrediction();
}

//____________________________________________________________________72
void NuCrossSectionFitter::SetupPars()
{

  //Set up parameters here
  Int_t parnum=0;
  for (UInt_t i=0; i<(fvNDTrueENuMuCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuCCBins.size()-1 == i){
      fFitter->SetParameter(parnum,
                          Form("NuMuOverflow",i),
                          1.0,0.2,1.0,0.0);
      fFitter->FixParameter(parnum);
    }
    else{
      fFitter->SetParameter(parnum,Form("NuMu%d",i),1.0,0.2,1.0,0.0);
    }
  }

  for (UInt_t i=0; i<(fvNDTrueENuMuBarCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuBarCCBins.size()-1 == i){
      fFitter->SetParameter(parnum,
                          Form("NuMuBarOverFlow",i),
                          1.0,0.2,1.0,0.0);
      fFitter->FixParameter(parnum);
    }
    else{
      fFitter->SetParameter(parnum,Form("NuMuBar%d",i),1.0,0.2,1.0,0.0);
    }
  }
  fFitter->SetParameter(parnum,"OtherPar",1.0,0.2,1.0,0.0);
  fFitter->FixParameter(parnum);
  parnum++;
  fFitter->SetParameter(parnum,"ShwEShift",1.0,0.0001,1.0,0.0);
  fFitter->FixParameter(parnum);
  parnum++;
  fFitter->SetParameter(parnum,"TrkEShift",1.0,0.2,1.0,0.0);
  fFitter->FixParameter(parnum);

  return;
}

//____________________________________________________________________72
void NuCrossSectionFitter::ReleaseNonOverflowNonNDFitPars()
{
  const TFitterMinuit* stupidRoot = fFitter;
  const vector<MinuitParameter>& vMinuitPars =
    stupidRoot->State().MinuitParameters();
  
  Int_t parnum=0;
  for (UInt_t i=0; i<(fvNDTrueENuMuCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuCCBins.size()-1 == i){
      if (!vMinuitPars[parnum].IsFixed()){
        fFitter->FixParameter(parnum);
      }
    }
    else{
      if (vMinuitPars[parnum].IsFixed()){
        fFitter->ReleaseParameter(parnum);
      }
    }
  }
  for (UInt_t i=0; i<(fvNDTrueENuMuBarCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuBarCCBins.size()-1 == i){
      if (!vMinuitPars[parnum].IsFixed()){
        fFitter->FixParameter(parnum);
      }
    }
    else{
      if (vMinuitPars[parnum].IsFixed()){
        fFitter->ReleaseParameter(parnum);
      }
    }
  }
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  parnum++;
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  parnum++;
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }

  return;
}

//____________________________________________________________________72
void NuCrossSectionFitter::ReleaseNonOverflowPars()
{
  const TFitterMinuit* stupidRoot = fFitter;
  const vector<MinuitParameter>& vMinuitPars =
    stupidRoot->State().MinuitParameters();
  
  Int_t parnum=0;
  for (UInt_t i=0; i<(fvNDTrueENuMuCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuCCBins.size()-1 == i){
      if (!vMinuitPars[parnum].IsFixed()){
        fFitter->FixParameter(parnum);
      }
    }
    else{
      if (vMinuitPars[parnum].IsFixed()){
        fFitter->ReleaseParameter(parnum);
      }
    }
  }
  for (UInt_t i=0; i<(fvNDTrueENuMuBarCCBins.size()); ++i,++parnum){
    if (fvNDTrueENuMuBarCCBins.size()-1 == i){
      if (!vMinuitPars[parnum].IsFixed()){
        fFitter->FixParameter(parnum);
      }
    }
    else{
      if (vMinuitPars[parnum].IsFixed()){
        fFitter->ReleaseParameter(parnum);
      }
    }
  }
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  parnum++;
  if (vMinuitPars[parnum].IsFixed()){
    fFitter->ReleaseParameter(parnum);
  }
  parnum++;
  if (vMinuitPars[parnum].IsFixed()){
    fFitter->ReleaseParameter(parnum);
  }

  return;
}

//____________________________________________________________________72
void NuCrossSectionFitter::FixAllParsExcept(const Int_t freeparnum)
{
  const TFitterMinuit* stupidRoot = fFitter;
  const vector<MinuitParameter>& vMinuitPars =
    stupidRoot->State().MinuitParameters();
  
  Int_t parnum=0;
  for (UInt_t i=0; i<(fvNDTrueENuMuCCBins.size()); ++i,++parnum){
    if (!vMinuitPars[parnum].IsFixed()){
      fFitter->FixParameter(parnum);
    }
  }
  for (UInt_t i=0; i<(fvNDTrueENuMuBarCCBins.size()); ++i,++parnum){
    if (!vMinuitPars[parnum].IsFixed()){
      fFitter->FixParameter(parnum);
    }
  }
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  parnum++;
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  parnum++;
  if (!vMinuitPars[parnum].IsFixed()){
    fFitter->FixParameter(parnum);
  }
  
  fFitter->ReleaseParameter(freeparnum);

  return;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::DoNDFit()
{  
  this->SetupPars();

  this->FixAllParsExcept(0);
  //Do the fit
  fFitter->Minimize();
  fFitter->PrintResults(0,0);
  
  this->FixAllParsExcept(1);
  //Do the fit
  fFitter->Minimize();
  fFitter->PrintResults(0,0);
  
  this->ReleaseNonOverflowNonNDFitPars();
  fFitter->Minimize();
  fFitter->PrintResults(0,0);
  /*
  this->ReleaseNonOverflowPars();
  fFitter->Minimize();
  fFitter->PrintResults(0,0);
  */
  const vector<Double_t> vBestFit = this->GetBestFitPars();
  this->FillTruthWeightedHistograms(vBestFit);
  fOutputFile->cd();
  fhNDMCRecoENuMuCC->Write("BestFitNDMCRecoENuMuCC");
  fhNDMCRecoENuMuBarCC->Write("BestFitNDMCRecoENuMuBarCC");
  return true;
}

//____________________________________________________________________72
const vector<Double_t> NuCrossSectionFitter::GetBestFitPars() const
{  
  vector<Double_t> vbestfit;
  Int_t parnum=0;
  for (UInt_t i=0; i<(fvNDTrueENuMuCCBins.size()); ++i,++parnum){
    vbestfit.push_back(fFitter->GetParameter(parnum));
  }
  for (UInt_t i=0; i<(fvNDTrueENuMuBarCCBins.size()); ++i,++parnum){
    vbestfit.push_back(fFitter->GetParameter(parnum));
  }
  vbestfit.push_back(fFitter->GetParameter(parnum));
  parnum++;
  vbestfit.push_back(fFitter->GetParameter(parnum));
  parnum++;
  vbestfit.push_back(fFitter->GetParameter(parnum));

  return vbestfit;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::Config()
{
  Bool_t goodConfig = true;
  if (!this->ConfigND()) goodConfig = false;
  if (!this->ConfigFD()) goodConfig = false;
  return goodConfig;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::ConfigND()
{
  Bool_t goodConfig = true;
  if (!this->CreateNDHistograms()) goodConfig = false;

  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "NDHistograms created"
    << endl;

  if (!this->FillNDDataHistograms()) goodConfig = false;

  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "NDHistograms filled"
    << endl;

  if (!this->CacheNDMCEvents()) goodConfig = false;

  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "ND MC events cached"
    << endl;

  return goodConfig;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::ConfigFD()
{
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDHistograms()
{
  Bool_t goodConfig = true;
  if (!this->CreateNDNuMuCCDataHisto()) goodConfig = false;
  if (!this->CreateNDNuMuBarCCDataHisto()) goodConfig = false;
  if (!this->CreateNDNuMuCCMCHisto()) goodConfig = false;
  if (!this->CreateNDNuMuBarCCMCHisto()) goodConfig = false;
  if (!this->CreateNDRawMCRecoECCNuMuHisto()) goodConfig = false;
  if (!this->CreateNDRawMCRecoECCNuMuBarHisto()) goodConfig = false;
  if (!this->CreateNDRawMCTrueECCNuMuHisto()) goodConfig = false;
  if (!this->CreateNDRawMCTrueECCNuMuBarHisto()) goodConfig = false;
  return goodConfig;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDNuMuCCMCHisto()
{
  if (fhNDMCRecoENuMuCC) {
    delete fhNDMCRecoENuMuCC; fhNDMCRecoENuMuCC = 0;
  }

  if (fhNDNuMuCCDataFromFile && fhNDNuMuBarCCDataFromFile){
    fhNDMCRecoENuMuCC = new TH1D(*fhNDNuMuCCDataFromFile);
    fhNDMCRecoENuMuCC->Reset();
    return true;
  }

  Int_t size = fvNDRecoENuMuCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuCCBins.begin();
       it != fvNDRecoENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhNDMCRecoENuMuCC = new TH1D("fhNDMCRecoENuMuCC",
                               "",
                               size-1,
                               recoEBinEdges);

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

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDRawMCRecoECCNuMuHisto()
{
  if (fhNDRawMCRecoECCNuMu) {
    delete fhNDRawMCRecoECCNuMu; fhNDRawMCRecoECCNuMu = 0;
  }

  if (fhNDNuMuCCDataFromFile && fhNDNuMuBarCCDataFromFile){
    fhNDRawMCRecoECCNuMu = new TH1D(*fhNDNuMuCCDataFromFile);
    fhNDRawMCRecoECCNuMu->Reset();
    return true;
  }

  Int_t size = fvNDRecoENuMuCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuCCBins.begin();
       it != fvNDRecoENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDRawMCRecoECCNuMu = new TH1D("fhNDRawMCRecoECCNuMu",
                                  "",
                                  size-1,
                                  recoEBinEdges);
  
  if (!fhNDRawMCRecoECCNuMu) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDRawMCRecoECCNuMuBarHisto()
{
  if (fhNDRawMCRecoECCNuMuBar) {
    delete fhNDRawMCRecoECCNuMuBar; fhNDRawMCRecoECCNuMuBar = 0;
  }

  if (fhNDNuMuCCDataFromFile && fhNDNuMuBarCCDataFromFile){
    fhNDRawMCRecoECCNuMuBar = new TH1D(*fhNDNuMuBarCCDataFromFile);
    fhNDRawMCRecoECCNuMuBar->Reset();
    return true;
  }

  Int_t size = fvNDRecoENuMuBarCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuBarCCBins.begin();
       it != fvNDRecoENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDRawMCRecoECCNuMuBar = new TH1D("fhNDRawMCRecoECCNuMuBar",
                                  "",
                                  size-1,
                                  recoEBinEdges);
  
  if (!fhNDRawMCRecoECCNuMuBar) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDRawMCTrueECCNuMuHisto()
{
  if (fhNDRawMCTrueECCNuMu) {
    delete fhNDRawMCTrueECCNuMu; fhNDRawMCTrueECCNuMu = 0;
  }

  Int_t size = fvNDTrueENuMuCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND true 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 =
         fvNDTrueENuMuCCBins.begin();
       it != fvNDTrueENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDRawMCTrueECCNuMu = new TH1D("fhNDRawMCTrueECCNuMu",
                                  "",
                                  size-1,
                                  recoEBinEdges);
  
  if (!fhNDRawMCTrueECCNuMu) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDRawMCTrueECCNuMuBarHisto()
{
  if (fhNDRawMCTrueECCNuMuBar) {
    delete fhNDRawMCTrueECCNuMuBar; fhNDRawMCTrueECCNuMuBar = 0;
  }

  Int_t size = fvNDTrueENuMuBarCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND true 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 =
         fvNDTrueENuMuBarCCBins.begin();
       it != fvNDTrueENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDRawMCTrueECCNuMuBar = new TH1D("fhNDRawMCTrueECCNuMuBar",
                                  "",
                                  size-1,
                                  recoEBinEdges);
  
  if (!fhNDRawMCTrueECCNuMuBar) return false;
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDNuMuBarCCMCHisto()
{
  if (fhNDMCRecoENuMuBarCC){
    delete fhNDMCRecoENuMuBarCC; fhNDMCRecoENuMuBarCC = 0;
  }

  if (fhNDNuMuCCDataFromFile && fhNDNuMuBarCCDataFromFile){
    fhNDMCRecoENuMuBarCC = new TH1D(*fhNDNuMuBarCCDataFromFile);
    fhNDMCRecoENuMuBarCC->Reset();
    return true;
  }

  Int_t size = fvNDRecoENuMuBarCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuBarCCBins.begin();
       it != fvNDRecoENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDMCRecoENuMuBarCC = new TH1D("fhNDMCRecoENuMuBarCC",
                                  "",
                                  size-1,
                                  recoEBinEdges);

  if (!fhNDMCRecoENuMuBarCC){return false;}
  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDNuMuCCDataHisto()
{
  if (fhNDDataRecoENuMuCC) {
    delete fhNDDataRecoENuMuCC; fhNDDataRecoENuMuCC = 0;
  }

  Int_t size = fvNDRecoENuMuCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuCCBins.begin();
       it != fvNDRecoENuMuCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhNDDataRecoENuMuCC = new TH1D("fhNDDataRecoENuMuCC",
                               "",
                                 size-1,
                               recoEBinEdges);

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

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CreateNDNuMuBarCCDataHisto()
{
  if (fhNDDataRecoENuMuBarCC){
    delete fhNDDataRecoENuMuBarCC; fhNDDataRecoENuMuBarCC = 0;
  }

  Int_t size = fvNDRecoENuMuBarCCBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND 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 =
         fvNDRecoENuMuBarCCBins.begin();
       it != fvNDRecoENuMuBarCCBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }
  
  fhNDDataRecoENuMuBarCC = new TH1D("fhNDDataRecoENuMuBarCC",
                                  "",
                                  size-1,
                                  recoEBinEdges);

  if (!fhNDDataRecoENuMuBarCC){return false;}
  return true;
}

//____________________________________________________________________72
void NuCrossSectionFitter::NDData(const string ndDataFileName)
{
  if (ftNDDataTree){delete ftNDDataTree; ftNDDataTree = 0;}
  ftNDDataTree = new NuTreeWrapper(ndDataFileName);

  //Get the total PoT
  fndDataPoT = ftNDDataTree->GetPoT();
  //fndDataPoT = 1.0;
  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "ndDataPoT: " << fndDataPoT
    << endl;
}

//____________________________________________________________________72
void NuCrossSectionFitter::NDNuMuBarCCData(const TH1D& ndData)
{
  if (fhNDNuMuBarCCDataFromFile){
    delete fhNDNuMuBarCCDataFromFile; fhNDNuMuBarCCDataFromFile = 0;
  }
  fhNDNuMuBarCCDataFromFile = new TH1D(ndData);
  fhNDNuMuBarCCDataFromFile->SetDirectory(0);
}

//____________________________________________________________________72
void NuCrossSectionFitter::NDNuMuCCData(const TH1D& ndData)
{
  if (fhNDNuMuCCDataFromFile){
    delete fhNDNuMuCCDataFromFile; fhNDNuMuCCDataFromFile = 0;
  }
  fhNDNuMuCCDataFromFile = new TH1D(ndData);
  fhNDNuMuCCDataFromFile->SetDirectory(0);
}

//____________________________________________________________________72
void NuCrossSectionFitter::NDMC(const string ndMCFileName)
{
  if (ftNDMCTree){delete ftNDMCTree; ftNDMCTree = 0;}
  ftNDMCTree = new NuTreeWrapper(ndMCFileName);

  //Get the total PoT
  fndMCPoT = ftNDMCTree->GetPoT();
  
  //fndMCPoT = 1.0;
  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "ndMCPoT: " << fndMCPoT
    << endl;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::CacheNDMCEvents()
{
  if (!ftNDMCTree){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND MC not supplied"
      << endl;
    return false;
  }

  fvNDMCCCNuMuEvents.clear();
  fvNDMCCCNuMuBarEvents.clear();

  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "Looping over ND MC:"
    << endl;
  Int_t mcEntries = ftNDMCTree->GetEntries();
  for (Int_t counter = 0;
       counter < mcEntries;
       ++counter){
    if (!(counter%10000)){
      MSG("NuCrossSectionFitter",Msg::kInfo)
        << counter << " of " << mcEntries
        << endl;
    }
    NuEvent currentEvent = ftNDMCTree->GetInfoObject(counter);
    fnuSystematic->Shift(currentEvent);
    NuSelectionResult_t recoType = this->SelectedAs(currentEvent);
    Double_t eventWeight = currentEvent.rw;
    if (NuXSecFit::kCCNuMu==recoType){
      fvNDMCCCNuMuEvents.push_back(currentEvent);
      fhNDRawMCRecoECCNuMu->Fill(currentEvent.energy,eventWeight);
      fhNDRawMCTrueECCNuMu->Fill(currentEvent.neuEnMC,eventWeight);
    }
    if (NuXSecFit::kCCNuMuBar==recoType){
      fvNDMCCCNuMuBarEvents.push_back(currentEvent);
      fhNDRawMCRecoECCNuMuBar->Fill(currentEvent.energy,eventWeight);
      fhNDRawMCTrueECCNuMuBar->Fill(currentEvent.neuEnMC,eventWeight);
    }
  }

  if (fndMCPoT){
    fhNDRawMCRecoECCNuMu->Scale(fndDataPoT/fndMCPoT);
    fhNDRawMCRecoECCNuMuBar->Scale(fndDataPoT/fndMCPoT);
    fhNDRawMCTrueECCNuMu->Scale(fndDataPoT/fndMCPoT);
    fhNDRawMCTrueECCNuMuBar->Scale(fndDataPoT/fndMCPoT);
  }

  fOutputFile->cd();
  fhNDRawMCRecoECCNuMu->Write("NDRawMCRecoECCNuMu");
  fhNDRawMCTrueECCNuMu->Write("NDRawMCTrueECCNuMu");
  fhNDRawMCRecoECCNuMuBar->Write("NDRawMCRecoECCNuMuBar");
  fhNDRawMCTrueECCNuMuBar->Write("NDRawMCTrueECCNuMuBar");

  return true;
}

//____________________________________________________________________72
const Bool_t NuCrossSectionFitter::FillNDDataHistograms()
{
  //Preferentially take data from histograms if they've been supplied
  if (fhNDNuMuCCDataFromFile && fhNDNuMuBarCCDataFromFile){
    if (fhNDDataRecoENuMuCC){
      delete fhNDDataRecoENuMuCC; fhNDDataRecoENuMuCC = 0;
    }
    if (fhNDDataRecoENuMuBarCC){
      delete fhNDDataRecoENuMuBarCC; fhNDDataRecoENuMuBarCC = 0;;
    }
    fhNDDataRecoENuMuCC = new TH1D(*fhNDNuMuCCDataFromFile);
    fhNDDataRecoENuMuBarCC = new TH1D(*fhNDNuMuBarCCDataFromFile);

    fOutputFile->cd();
    fhNDDataRecoENuMuCC->Write("InputNDNuMuCCData");
    fhNDDataRecoENuMuBarCC->Write("InputNDNuMuBarCCData");
    return true;
  }

  if (!fhNDDataRecoENuMuCC){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND numu data histogram not created"
      << endl;
    return false;
  }
  if (!fhNDDataRecoENuMuBarCC){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND numubar data histogram not created"
      << endl;
    return false;
  }
  if (!ftNDDataTree){
    MSG("NuCrossSectionFitter",Msg::kWarning)
      << "ND data not supplied"
      << endl;
    return false;
  }

  fhNDDataRecoENuMuCC->Reset();
  fhNDDataRecoENuMuBarCC->Reset();

  MSG("NuCrossSectionFitter",Msg::kInfo)
    << "Looping over ND data:"
    << endl;
  Int_t dataEntries = ftNDDataTree->GetEntries();
  for (Int_t counter = 0;
       counter < dataEntries;
       ++counter){
    if (!(counter%10000)){
      MSG("NuCrossSectionFitter",Msg::kInfo)
        << counter << " of " << dataEntries
        << endl;
    }
    NuEvent currentEvent = ftNDDataTree->GetInfoObject(counter);
    NuSelectionResult_t recoType = this->SelectedAs(currentEvent);
    Double_t recoEnergy = currentEvent.energy;
    
    Int_t truthpar = this->FakeTruthIndex(currentEvent);
    Double_t nuweight = 1.0;
    nuweight *= currentEvent.rw;
    if (truthpar>=0){
      nuweight *= fvfakepars[truthpar];
    }

    if (recoType==NuXSecFit::kCCNuMu){
      fhNDDataRecoENuMuCC->Fill(recoEnergy,nuweight);
    }
    if (recoType==NuXSecFit::kCCNuMuBar){
      fhNDDataRecoENuMuBarCC->Fill(recoEnergy,nuweight);
    }
  }

  fOutputFile->cd();
  fhNDDataRecoENuMuCC->Write("InputNDNuMuCCData");
  fhNDDataRecoENuMuBarCC->Write("InputNDNuMuBarCCData");

  return true;
}

//____________________________________________________________________72
const NuSelectionResult_t
NuCrossSectionFitter::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;
}

---