NCExtrapolationRS.cxx

Go to the documentation of this file.

// $Id: NCExtrapolationRS.cxx,v 1.31 2008/01/11 15:53:31 bckhouse Exp $
//
// A class to implement Tobi's x-section NearFit and use the results to
// produce FD predictions near to far extrapolations
//
// A. Sousa 5/2007

#include "NCUtils/Extrapolation/NCExtrapolationRS.h"
#include "NCUtils/Extrapolation/NCEnergyBin.h"
#include "MessageService/MsgService.h"
#include "MCReweight/ReweightHelpers.h"
#include "AnalysisNtuples/ANtpAnalysisInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"

#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TLegend.h"
#include "TLine.h"

#include <cassert>

CVSID("$Id: NCExtrapolationRS.cxx,v 1.31 2008/01/11 15:53:31 bckhouse Exp $");

NCExtrapolationRS* NCExtrapolationRS::sFit = NULL;

//......................................................................
NCExtrapolationRS::NCExtrapolationRS() :
  NCExtrapolation::NCExtrapolation()
{

}

//.......................................................................
NCExtrapolationRS::NCExtrapolationRS(NCAnalysisUtils *utils,
                                     std::map<Int_t,double>& dataNear,
                                     std::map<Int_t,double>& mcNear,
                                     std::map<Int_t,double>& dataFar,
                                     std::map<Int_t,double>& mcFar,
                                     bool useCC,
                                     bool useNC) :
  NCExtrapolation(utils, dataNear, mcNear, dataFar, mcFar, useCC, useNC),
  fUseChi2Function(false),
  fAnnCutMC(0.675),
  fAnnCutData(0.675),
  fCutLowBins(false),
  fMoveMinStrip(0),
  fRandom(0)
{
  MSG("NCExtrapolationRS",Msg::kDebug) << "NCExtrapolationRS constructor"
                                       << endl;
  fExtrapolationType = NCType::kExtrapolationRS;

  //fill the map of data to pot ratios for the near detector
  std::map<Int_t, Double_t>::const_iterator potItr = dataNear.begin();
  fDataFraction.clear();
  while(potItr != dataNear.end()){
    fDataFraction[potItr->first] = potItr->second;
    ++potItr;
  }

  potItr = mcNear.begin();
  while(potItr != mcNear.end()){
    if(potItr->second == 0.)
      fDataFraction[potItr->first] *= 0.;
    else
      fDataFraction[potItr->first] /= potItr->second;

    ++potItr;
  }

  // book all the necessary histograms
  this->BookHistograms();

  //fBeams is not yet created, so use scaleFactor maps instead.
  for (UInt_t i=0; i<fDataFraction.size(); i++) {  // loop over different
                                                      // beams
    // create vectors to hold MC and Data events
    fMCListNear.push_back  ( new std::vector<MCEvent> );
    fDataListNear.push_back( new std::vector<MCEvent> );
    fMCListFar.push_back   ( new std::vector<MCEvent> );
    fDataListFar.push_back ( new std::vector<MCEvent> );

    fScaleToData.push_back(true);

  } // end beams loop

  // this is needed for fit and static fit function.
  sFit = this;
}

//----------------------------------------------------------------------
NCExtrapolationRS::~NCExtrapolationRS()
{
}

//----------------------------------------------------------------------
//this method allows the user to pass in appropriate settings for the
//extrapolation such as locations of files, whether to generate the files
//again, etc.  the registry comes from the job module GetConfig method
void NCExtrapolationRS::Prepare(const Registry& r)
{
  MSG("NCExtrapolationRS", Msg::kWarning)
    << "NCExtrapolationRS::Prepare(const Registry& r)" << endl;
//   int         tmpb;  // a temp bool. See comment under DefaultConfig...
//   int         tmpi;  // a temp int.
//   double      tmpd;  // a temp double.
//   const char* tmps;  // a temp string

  NCExtrapolation::Prepare(r);

  int tmpb;
  const char* tmps;  // a temp string

  // check if we are running on a limited number of test events
  // --> don't scale with total POT in this case
  r.Get("CountChainEntries",tmpb);
  if (!tmpb) {
    for (uint i=0; i< fScaleToData.size(); ++i)
      fScaleToData.at(i) = false;
  }

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

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationRS::AddEvent(ANtpHeaderInfo *headerInfo,
                                 ANtpBeamInfo *beamInfo,
                                 ANtpRecoInfo *recoInfo,
                                 ANtpAnalysisInfo *analysisInfo,
                                 ANtpTruthInfoBeam *truthInfo,
                                 int runType,
                                 bool useMCAsData,
                                 int fileType)
{

  if(headerInfo->detector == (int)Detector::kNear)
    AddNearEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                 truthInfo, runType, useMCAsData, fileType);

  if(headerInfo->detector == (int)Detector::kFar)
    AddFarEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                truthInfo, runType, useMCAsData, fileType);

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationRS::AddNearEvent(ANtpHeaderInfo *headerInfo,
                                     ANtpBeamInfo *beamInfo,
                                     ANtpRecoInfo *recoInfo,
                                     ANtpAnalysisInfo *analysisInfo,
                                     ANtpTruthInfoBeam *truthInfo,
                                     int runType,
                                     bool useMCAsData,
                                     int fileType)
{
  BeamType::BeamType_t beamType = BeamType::FromZarko(beamInfo->beamType);
  NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                            truthInfo, runType, useMCAsData, fileType);
  int beam = NCExtrapolation::FindBeamIndex(Detector::kNear, beamType, runType);

  if (HistEvisNCd.at(beam)==0) {
    MSG("NCExtrapolationRS",Msg::kError) << "Histograms with index "
                                         << beam << ", Beam "
                                         << BeamType::AsString(beamType)
                                         << ", don't exist. Exiting..."
                                         << endl;
    exit(1);
  }else{
    MCEvent event;

//  event.Evis   = recoInfo->nuEnergyNC; // reconstructed energy
//  event.EvisMc = recoInfo->nuEnergyNC;
    event.Evis   = recoInfo->nuEnergy; // reconstructed energy
    event.EvisMc = recoInfo->nuEnergy;

    event.EvisCC = recoInfo->nuEnergy;   // CC energy incl. track
    event.EshowerCC = recoInfo->showerEnergyCC;
    event.EshowerNC = recoInfo->showerEnergyNC;
    event.Etrack = recoInfo->muEnergy;   // track energy
    if (event.Etrack<0) event.Etrack = 0;

    event.Weight = 1.;          // this is a systematic weight
    // for the event

    if (truthInfo) event.BeamWeight=recoInfo->weight;
    else event.BeamWeight=1.0;

    //Handle MC event
    if(truthInfo && !useMCAsData) {
      event.Enu    = truthInfo->nuEnergy;    // neutrino energy
      event.McType = truthInfo->interactionType;
      event.nuDCosX = truthInfo->nuDCosX;
      event.nuDCosY = truthInfo->nuDCosY;
      event.nuDCosZ = truthInfo->nuDCosZ;
      event.targetEnergy = truthInfo->targetEnergy;
      event.targetPX = truthInfo->targetPX;
      event.targetPY = truthInfo->targetPY;
      event.targetPZ = truthInfo->targetPZ;
      event.nuFlavor = truthInfo->nuFlavor;
      event.interactionType = truthInfo->interactionType;
      event.iresonance = truthInfo->resonanceCode;
      event.hadronicY = truthInfo->hadronicY;
      event.X = truthInfo->bjorkenX;
      event.q2 = truthInfo->q2;
      event.w2 = truthInfo->w2;

      Int_t nucleus = ReweightHelpers::FindNucleusNumber(
        static_cast<Int_t>(truthInfo->atomicNumber),
        static_cast<Int_t>(truthInfo->atomicWeight)
        );
      event.nucleus = nucleus;
      event.mcstate = truthInfo->initialState;
      event.hadronicFinalState = truthInfo->hadronicFinalState;

    }//if(truthInfo)

    // it is not due to oscillation
    //    event.SelType = analysisTRann->isCC;
    // by moving fANNCutMC, we investigate systematics
    //if (analysisTRann->separationParameter < fAnnCutData)
    if(analysisInfo->isNC > 0)
      event.SelType = 0;      // this is NC
    else if (analysisInfo->isCC > 0)
      event.SelType = 1;      // this is CC
    else
      event.SelType = -1;
    // cut out non-fid and un-clean events  --> selType = -1
    if (recoInfo->inFiducialVolume==0 || recoInfo->isCleanHighMultSnarl==0)
      event.SelType = -1;

    // now fill histograms with real data or MC

    if (truthInfo && !useMCAsData){

      if (event.SelType!=-1) {         //fiducial and cleaning
        HistEnu.at(beam)->Fill(event.Enu);

        if (event.SelType==1)  {            // selected as CC
          HistEvisCC.at(beam)->Fill(event.EvisCC);
          HistEvisCCwt.at(beam)->Fill(event.EvisCC,event.BeamWeight);

          if (event.McType==0) {
            HistEvisCCb.at(beam)->Fill(event.EvisCC, event.BeamWeight);
          }
        }
        else if (event.SelType==0){         // selected as NC
          HistEvisNC.at(beam)->Fill(event.Evis);
          HistEvisNCwt.at(beam)->Fill(event.Evis,event.BeamWeight);

          if (event.McType==1) {
            HistEvisNCb.at(beam)->Fill(event.Evis, event.BeamWeight);
          }
        }
        fMCListNear.at(beam)->push_back(event);
      }
    }//if(truthInfo && !useMCAsData)
    else{ //data event

      if (event.SelType!=-1) {
        if (event.SelType==0) {
          MAXMSG("NCExtrapolationRS",Msg::kInfo,5) << "BeamWeight: "
                                                   << event.BeamWeight
                                                   << endl;
          HistEvisNCd.at(beam)->Fill(event.Evis, event.BeamWeight);
        }
        if (event.SelType==1){
          HistEvisCCd.at(beam)->Fill(event.EvisCC, event.BeamWeight);
        }
      }

      fDataListNear.at(beam)->push_back(event);

    }//else(truthInfo)

  } //else(HistEvis)


  return;
}
//----------------------------------------------------------------------
void NCExtrapolationRS::AddFarEvent(ANtpHeaderInfo *headerInfo,
                                    ANtpBeamInfo *beamInfo,
                                    ANtpRecoInfo *recoInfo,
                                    ANtpAnalysisInfo *analysisInfo,
                                    ANtpTruthInfoBeam *truthInfo,
                                    int runType,
                                    bool useMCAsData,
                                    int fileType)
{
  MAXMSG("NCExtrapolationRS",Msg::kInfo,1) << "Running AddFarEvent()..."
                                           << endl;

  // find correct trueBin to apply fit result scale factor;
  Int_t trueBin = -1;
  if (truthInfo) {
    Float_t Enu = truthInfo->nuEnergy;
    if (Enu < 4)
      trueBin = 0;
    else if (Enu < 8)
      trueBin = 1;
    else if (Enu < 15)
      trueBin = 2;
    else if (Enu < 120)
      trueBin = 3;
    else
      trueBin = 99;
  }

  if(truthInfo && !useMCAsData && fileType != NCType::kTauFile) {
    if (analysisInfo->isNC > 0) {
      fNCVecFar.push_back(recoInfo->showerEnergy);

//       now filled in NCBeam::AddEvent() method
//       fBeams[beam].GetDefaultMCHistogram(NCType::kNC)->Fill(recoInfo->nuEnergy,
//                                                          recoInfo->weight);

      // scale FD monte carlo by ND fit results
      if (truthInfo->interactionType==0 && trueBin!=99) {
        recoInfo->weight *= mScaleNC.at(trueBin);

        if (mScaleNC.at(trueBin)!= 1)
          cout << "Scale factor: " << mScaleNC.at(trueBin) << endl;
      }
      fNCWeightFar.push_back(recoInfo->weight);
    }
    else if(analysisInfo->isCC > 0) {
      fCCVecFar.push_back(recoInfo->showerEnergy + recoInfo->muEnergy);

//       now filled in NCBeam::AddEvent() method
//       fBeams[beam].GetDefaultMCHistogram(NCType::kCC)->Fill(recoInfo->nuEnergy,
//                                                          recoInfo->weight);

      // scale FD monte carlo by ND fit results
      Double_t weight = recoInfo->weight;
      if (truthInfo->interactionType==0 && trueBin!=99)
        recoInfo->weight *= mScaleNC.at(trueBin);
      fCCWeightFar.push_back(weight);

    }
  }

  NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
                            truthInfo, runType, useMCAsData, fileType);

  return;
}

//----------------------------------------------------------------------
//Book histograms to be used in fit.
void NCExtrapolationRS::BookHistograms()
{

  TFile *nu = new TFile("$SRT_LOCAL/NCUtils/data/FluxErrorHistos14.root",
                        "READ");
  TFile *nubar = new TFile("$SRT_LOCAL/NCUtils/data/FluxErrorHistos-14.root",
                           "READ");

  if(!nu->IsOpen() || !nubar->IsOpen()){
    MAXMSG("NCExtrapolationRS",Msg::kInfo, 20)<< "SKZP error files "
                                           << "not found! "
                                           << "Not filling SKZP histos..."
                                           << endl;
  }

  std::vector<TH1D*> HistCCArray;
  std::vector<TH1D*> HistNCArray;

  MSG("NCExtrapolationRS",Msg::kDebug) << "fDataFraction.size(): "
                                       << fDataFraction.size()
                                       << endl;

  gROOT->cd();
  for(uint i = 0; i < fDataFraction.size(); ++i){
    string name = (string) Form("HistEnu_%i",i);
    HistEnu.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisCC_%i", i);
    HistEvisCC.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisNC_%i", i);
    HistEvisNC.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisCCwt_%i", i);
    HistEvisCCwt.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisNCwt_%i", i);
    HistEvisNCwt.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisCCrew_%i", i);
    HistEvisCCrew.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisNCrew_%i", i);
    HistEvisNCrew.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisCCb_%i", i);
    HistEvisCCb.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisNCb_%i", i);
    HistEvisNCb.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisCCd_%i", i);
    HistEvisCCd.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    name = (string) Form("HistEvisNCd_%i", i);
    HistEvisNCd.push_back( new TH1D(name.c_str(), name.c_str(), 60, 0., 30.) );

    for (Int_t trueBin = 0; trueBin<4; trueBin++) {

      name = (string) Form("HistEvisCCtrueBin_%i_%i", i, trueBin);
      HistCCArray.push_back( new TH1D(name.c_str(), name.c_str(),
                                      60, 0., 30.) );

      name = (string) Form("HistEvisNCtrueBin_%i_%i", i, trueBin);
      HistNCArray.push_back( new TH1D(name.c_str(), name.c_str(),
                                      60, 0., 30.) );

      // also set scale factors to 1 here -- this is the default
      mInputScale.push_back(1.0);
    }

    HistEvisCCtrueBin.push_back(HistCCArray);
    HistEvisNCtrueBin.push_back(HistNCArray);

    HistCCArray.clear();
    HistNCArray.clear();

    // SKZP flux error histograms
    name = (string) Form("HistCCFluxErrorMinus_%i", i);
    HistCCFluxErrorMinus.push_back( new TH1D(name.c_str(), name.c_str(),
                                             60, 0., 30.) );
    name = (string) Form("HistNCFluxErrorMinus_%i", i);
    HistNCFluxErrorMinus.push_back( new TH1D(name.c_str(), name.c_str(),
                                             60, 0., 30.) );
    name = (string) Form("HistCCFluxErrorPlus_%i", i);
    HistCCFluxErrorPlus.push_back( new TH1D(name.c_str(), name.c_str(),
                                            60, 0., 30.) );
    name = (string) Form("HistNCFluxErrorPlus_%i", i);
    HistNCFluxErrorPlus.push_back( new TH1D(name.c_str(), name.c_str(),
                                            60, 0., 30.) );

    if(nu->IsOpen() && nubar->IsOpen()){
      //read SKZP flux error histograms from file
      name = (string) Form("HistTrueEFluxErrPlus_nu_%i", i);
      HistTrueEFluxErrPlus_nu.push_back( dynamic_cast<TH1D*>
                                         ( nu->Get(name.c_str() )) );
      if (! HistTrueEFluxErrPlus_nu[i]) {
        MSG("NExtrapolationRS",Msg::kError) << "Histogram "
                                            << name
                                            << " not found!"
                                            << endl;
        exit(1);
      }
      name = (string) Form("HistTrueEFluxErrMinus_nu_%i", i);
      HistTrueEFluxErrMinus_nu.push_back( dynamic_cast<TH1D*>
                                          ( nu->Get(name.c_str())) );
      if (! HistTrueEFluxErrMinus_nu[i]) {
        MSG("NCextrapolationRS",Msg::kError) << "Histogram "
                                             << name
                                             << " not found!"
                                             << endl;
        exit(1);
      }

      name = (string) Form("HistTrueEFluxErrPlus_nubar_%i", i);
      HistTrueEFluxErrPlus_nubar.push_back( dynamic_cast<TH1D*>
                                            ( nubar->Get(name.c_str())) );
      if (! HistTrueEFluxErrPlus_nubar[i]) {
        MSG("NCExtrapolationRS",Msg::kError) << "Histogram "
                                             << name
                                             << " not found!"
                                             << endl;
        exit(1);
      }

      name = (string) Form("HistTrueEFluxErrMinus_nubar_%i", i);
      HistTrueEFluxErrMinus_nubar.push_back( dynamic_cast<TH1D*>
                                             (nubar->Get(name.c_str())) );
      if (! HistTrueEFluxErrMinus_nubar[i]) {
        MSG("NCExtrapolationRS",Msg::kError) << "Histogram "
                                             << name
                                             << " not found!"
                                             << endl;
        exit(1);
      }
    }//if skzp file open

  } // fBeams.size();

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationRS::DoFit(Int_t PrintLevel)
{
//  PlotData();

  // give start values and errors for fit
  // parameters to scale the NC xsec -- use 4 bins
  // initialize the vector for final numbers "mScaleNC" and "mScaleNCErr"
  // with starting values as well
  for (Int_t i=0; i<4; i++) {
    m0ScaleNC.push_back(1.0);
//    m0ScaleNCErr.push_back(0.005);
    m0ScaleNCErr.push_back(0.);
    mScaleNC.push_back(1.0);
//    mScaleNCErr.push_back(0.005);
    mScaleNCErr.push_back(0.);
  }
  // and now do the fit
  Fit(PrintLevel);

  // result of fit is now in mXxx and mXxxErr

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationRS::Fit(Int_t PrintLevel)
{

  // set up a Minuit to simultanisouly fit CC & NC spectra

  Double_t arglist[10];
  Int_t    ierr;

  TMinuit* gMinuit = new TMinuit(4);

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

  arglist[0] = PrintLevel;
  gMinuit->mnexcm("SET PRInt",arglist,1,ierr);
  if ( PrintLevel<-1 )
    gMinuit->mnexcm("SET NOWarning",arglist,0,ierr);
//  arglist.at(0) = 0.5;
  arglist[0] = 1; // added a factor 2 to the likelihood function
  gMinuit->mnexcm("SET ERRor",arglist,1,ierr);
  //arglist.at(0) = 1e-5; // precision of function
  //gMinuit->mnexcm("SET EPS",arglist,1,ierr);
  arglist[0] = 1; // minimization strategy
  gMinuit->mnexcm("SET STR",arglist,1,ierr);

  // Set starting values
  Double_t par[4];
  for (Int_t i=0; i<4; i++)
    par[i] = m0ScaleNC.at(i);

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

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

  gMinuit->FixParameter(0);
  gMinuit->FixParameter(1);
  gMinuit->FixParameter(2);

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

  gMinuit->Migrad();

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

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

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

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

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

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

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

  // plot result

  if ( PrintLevel >-2 ) {

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

    for (Int_t i=0; i<4; i++)
      pars[i] = mScaleNC.at(i);

    LogLikelihoodFunc(npar, pars, LL, pars, 0);


    std::vector<TH1D*> ratCC(fDataFraction.size(),0); // pointers for ratio
                                                      //  histogram
    std::vector<TH1D*> ratCCfit(fDataFraction.size(),0);
    std::vector<TH1D*> ratNC(fDataFraction.size(),0);
    std::vector<TH1D*> ratNCfit(fDataFraction.size(),0);

    for (UInt_t i=0; i<fDataFraction.size(); i++) {
      gROOT->cd();
      // Draw the best fit spectrum
      string canvas_name = (string) Form("cSpec_%i",i);
      TCanvas* cSpec = dynamic_cast<TCanvas*>
        ( gROOT->FindObjectAny(canvas_name.c_str()) );
      if (!cSpec) cSpec = new TCanvas(canvas_name.c_str(),canvas_name.c_str(),
                                      40,40,500,700);
      cSpec->Divide(1,2);
      cSpec->cd(1);
      HistEvisCCwt.at(i)->Draw();
      HistEvisCCrew.at(i)->SetLineColor(6);
      HistEvisCCrew.at(i)->Draw("sames");
      HistEvisCCd.at(i)->Draw("samepe");

      cSpec->cd(2);
      HistEvisNCwt.at(i)->Draw();
      HistEvisNCrew.at(i)->SetLineColor(6);
      HistEvisNCrew.at(i)->Draw("sames");
      HistEvisNCd.at(i)->Draw("samepe");

      // Draw the best fit mc/data ration
      canvas_name = (string) Form("cRat_%i",i);
      TCanvas* cRat = dynamic_cast<TCanvas*>
        ( gROOT->FindObjectAny(canvas_name.c_str()) );
      if (!cRat)  cRat = new TCanvas(canvas_name.c_str(),canvas_name.c_str(),
                                     40,40,500,700);
      cRat->Divide(1,2);
      ratCC.at(i) = dynamic_cast<TH1D*>( HistEvisCCd.at(i)->Clone("ratCC") );
      ratCCfit.at(i) = dynamic_cast<TH1D*>(
        HistEvisCCd.at(i)->Clone("ratCCfit") );
      ratCC.at(i)->Divide( HistEvisCCwt.at(i) );
      ratCCfit.at(i)->Divide( HistEvisCCrew.at(i) );
      ratNC.at(i) = dynamic_cast<TH1D*>( HistEvisNCd.at(i)->Clone("ratNC") );
      ratNCfit.at(i) = dynamic_cast<TH1D*>(
        HistEvisNCd.at(i)->Clone("ratNCfit") );
      ratNC.at(i)->Divide( HistEvisNCwt.at(i) );
      ratNCfit.at(i)->Divide( HistEvisNCrew.at(i) );

      cRat->cd(1);
      ratCC.at(i)->Draw();
      ratCCfit.at(i)->SetLineColor(6);
      ratCCfit.at(i)->GetYaxis()->SetRangeUser(0.5,1.5);
      ratCCfit.at(i)->Draw("sames");
      TLine *line = new TLine(0,1.0,30,1.0);
      line->Draw("same");
      cRat->cd(2);
      ratNC.at(i)->Draw();
      ratNCfit.at(i)->SetLineColor(6);
      ratNCfit.at(i)->GetYaxis()->SetRangeUser(0.5,1.5);
      ratNCfit.at(i)->Draw("sames");
      line->Draw("same");
    }

    MSG("NCExtrapolationRS",Msg::kInfo)
      << "The final likelihood is: LL=" << setprecision(10) << LL << endl;

    // write some stuff to a log file
    *logFile << "Fit results:" << endl;
    logFile->setf(ios::fixed);
    logFile->precision(5);
    for (Int_t i=0; i<4; i++)
      *logFile << "Par " << i << ": " << mScaleNC[i] << " +"
               << mScalePos[i] << " " << mScaleNeg[i] << " +/- "
               << mScaleNCErr.at(i) << ", Corr: " << corr[i] << endl;

    *logFile << "Convergence Status: " << convergence << endl
             << "Final Likelihood LL=" << LL << endl;
      // << "Written histograms to " << histofile->GetName() << endl;
    logFile->unsetf(ios::fixed);


    // get the number of d.o.f. from chi2 fit
    Int_t dof = 0;
    if (fUseChi2Function) {
      for (UInt_t i = 0; i<fDataFraction.size(); ++i ) {
        Double_t chi2 = -1;
        Int_t ndf =  -1;
        Int_t igood = 0;
        HistEvisNCd.at(i)->Chi2TestX(HistEvisNCrew.at(i),chi2,ndf,igood,
                                     "UUCHI2");
        dof += ndf;
        MSG("NCExtrapolationRS",Msg::kDebug) << "Ndf for beam " << i << ": "
                                             << ndf << endl;
        if (igood != 0)
          MSG("NCExtrapolationRS",Msg::kWarning) << "Chi2 igood parameter"
                                                 <<" is non-zero: "
                                       << igood << endl;
      }
    }


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

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

  return;
}
//----------------------------------------------------------------------
void NCExtrapolationRS::FillFluxRecoHists(Int_t i, MCEvent event) {

  // get SKZP error from relevant histograms
  Double_t errPlus, errMinus;
  Int_t maxBin = HistTrueEFluxErrPlus_nu.at(i)->GetNbinsX();
  Int_t EtrueBin = HistTrueEFluxErrPlus_nu.at(i)->FindBin(event.Enu);
  // energies higher than 20GeV get the value for 20 GeV
  EtrueBin = min(EtrueBin, maxBin);

  if (event.nuFlavor==14) {
    errPlus = HistTrueEFluxErrPlus_nu.at(i)->GetBinContent(EtrueBin);
    errMinus = HistTrueEFluxErrMinus_nu.at(i)->GetBinContent(EtrueBin);
  }
  else if (event.nuFlavor==-14) {
    errPlus = HistTrueEFluxErrPlus_nubar.at(i)->GetBinContent(EtrueBin);
    errMinus = HistTrueEFluxErrMinus_nubar.at(i)->GetBinContent(EtrueBin);
  }
  else
    errPlus = errMinus = 1.0;

  if (event.SelType==1) {
    HistCCFluxErrorPlus.at(i)->Fill(event.EvisCC,event.Weight*
                                    errPlus*event.BeamWeight);
    HistCCFluxErrorMinus.at(i)->Fill(event.EvisCC,event.Weight*
                                     errMinus*event.BeamWeight);
  }
  else if (event.SelType==0) {
    HistNCFluxErrorPlus.at(i)->Fill(event.Evis, event.Weight*
                                    errPlus*event.BeamWeight);
    HistNCFluxErrorMinus.at(i)->Fill(event.Evis, event.Weight*
                                     errMinus*event.BeamWeight);
  }

  return;
}
//----------------------------------------------------------------------

Double_t NCExtrapolationRS::GetChiSquare(Int_t selection, Int_t beam) {

  // simple ChiSquare test of agreement between MC and data

  Double_t chi2 = 0;
  Int_t dof = 0;

  if (selection==1) {
    dof  = sFit->HistEvisCCd.at(beam)->GetNbinsX();

    for (Int_t i=1; i<=(dof+1); i++) {   //include overflow
      Double_t data = sFit->HistEvisCCd.at(beam)->GetBinContent(i);
      Double_t mc = sFit->HistEvisCCrew.at(beam)->GetBinContent(i);
      Double_t diff = data - mc;
      if ((mc > 5) && (data>0)) {
        chi2 += diff*diff/data;
      }
      else {
        MSG("NCExtrapolationRS",Msg::kWarning) << "Too few entries in bin "
                                               << i << ", beam " << beam
                                               << endl;
        dof--;
      }
    }
  }

  else if (selection==0) {
    dof  = sFit->HistEvisNCd.at(beam)->GetNbinsX();

    for (Int_t i=1; i<=(dof+1); i++) { // include overflow
      Double_t data = sFit->HistEvisNCd.at(beam)->GetBinContent(i);
      Double_t mc = sFit->HistEvisNCrew.at(beam)->GetBinContent(i);
      Double_t diff = data - mc;
      if ((mc > 5) && (data>0)) {
        chi2 += diff*diff/data;
      }
      else {
        MSG("NCExtrapolationRS",Msg::kWarning) << "Too few entries in bin "
                                               << i << ", beam " << beam
                                               << endl;
        dof--;
      }
    }
  }
  else {
    MSG("NCExtrapolationRS",Msg::kInfo)
      << "Invalid selection type. Chose 0==NC or 1==CC." << endl;
    return 0;
  }

  MSG("NCExtrapolationRS",Msg::kDebug)
    << "Calculated Chi2 " << chi2 << " for " << dof << " d.o.f."
    << endl;

  return chi2;
}

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

//LogLikelihood fit function.
void NCExtrapolationRS::LogLikelihoodFunc(Int_t &npar,
                                          Double_t *gin,
                                          Double_t &f,
                                          Double_t *pars,
                                          Int_t iflag)
{
  iflag++; // NOT USED FOR ANYTHING, JUST TO QUIET THE COMPILER
  if (0)  cout << gin << npar << endl;

  f = 0;

//   //  cout << "Parameter 5: " << pars[4];

//   // cross section reweighting MC using neugen parameters
//   //  sFit->Reweight(ma_qe,ma_res,kno112,kno122,kno113,kno123,false);

  for(uint i = 0; i < (sFit->fDataFraction).size(); ++i){
    sFit->HistEvisNCrew.at(i)->Reset();
    sFit->HistEvisCCrew.at(i)->Reset();
  }

  // now fill Reweighted Histograms
  // do something faster than histogram filling

  for(uint i = 0; i < (sFit->fBeams).size(); ++i){

    // do not use FD events in NearFit -- FD are also in fBeams vector
    if ((sFit->fBeams).at(i).GetDetector()!= Detector::kNear)
      continue;

    Int_t    CCStartBin = 0;
    Int_t    CCStopBin  = sFit->HistEvisCC.at(i)->GetNbinsX()+1;
    Double_t CCEmin     = sFit->HistEvisCC.at(i)->GetBinLowEdge(1);
    Double_t CCDeltaE   = sFit->HistEvisCC.at(i)->GetBinWidth(1);
    std::vector<Double_t> CCSignal(CCStopBin+1,0.);
    std::vector<Double_t> CCBackground(CCStopBin+1,0.);

    Int_t    NCStartBin = 0;
    Int_t    NCStopBin  = sFit->HistEvisNC.at(i)->GetNbinsX()+1;
    Double_t NCEmin     = sFit->HistEvisNC.at(i)->GetBinLowEdge(1);
    Double_t NCDeltaE   = sFit->HistEvisNC.at(i)->GetBinWidth(1);
    std::vector<Double_t> NCSignal(NCStopBin+1,0.);
    std::vector<Double_t> NCBackground(NCStopBin+1,0.);

    Int_t index = 0;

    std::vector<MCEvent>::iterator mcIt;
    for (mcIt=sFit->fMCListNear.at(i)->begin();
         mcIt<sFit->fMCListNear.at(i)->end(); mcIt++) {

      MCEvent event = (*mcIt);
      // use four wide bins in true E
      Int_t trueBin = -1;
      if (event.Enu < 4)
        trueBin = 0;
      else if (event.Enu < 8)
        trueBin = 1;
      else if (event.Enu < 15)
        trueBin = 2;
      else if (event.Enu < 120)
        trueBin = 3;
      else
        trueBin = 99;

      // include SKZP reweighting
      Double_t weight = event.BeamWeight*event.Weight;
      //Double_t weight = 1.0;

      if ( event.SelType==0 ) {
        // this is NC selected
        Double_t evis = event.Evis; // total calorimetric energy
        if (evis>100) index = NCStopBin;
        else {
          index = static_cast<Int_t>( ( evis-NCEmin )/ NCDeltaE ) + 1 ;
          index = min(index,NCStopBin);
          index = max(index,NCStartBin);
        }
        // do the scaling on true energy here
        if (event.McType==0 && trueBin<4)
          weight *= pars[trueBin];

        NCSignal.at(index) += weight;

        if (event.McType==1 ) {
          NCBackground.at(index) += weight;
        }
      }
      else if (event.SelType==1){
        // this in CC selected
        Double_t evis = event.EvisCC; // track+shower energy
        if (evis>100) index = CCStopBin;
        else {
          index = (Int_t) ( ( evis-CCEmin )/ CCDeltaE ) + 1;
          index = min(index,CCStopBin);
          index = max(index,CCStartBin);
        }
        if (event.McType==0 && trueBin<5) {
          // do the scaling on true energy here
          weight *= pars[trueBin];
        }
        CCSignal.at(index) += weight;

        if (event.McType ==0 ) {
          CCBackground.at(index) += weight;
        }
      }
    }

    Double_t CCInt     = 0.;
    Double_t NCInt     = 0.;
    Double_t CCdInt     = 0.;
    Double_t NCdInt     = 0.;
    Double_t binContent = 0.;
    //    Double_t scale = sFit->fDataFraction.at(i);
    Double_t scale =
      sFit->fDataFraction[BeamType::ToZarko((sFit->fBeams[i]).GetBeamType())];
    // decide whether data is scaled to MC or vice versa
    if (!sFit->fScaleToData.at(i))
      scale = 1.0;

    const Double_t minContent = .001;

    // fill histograms from arrays
    for (Int_t ibin = CCStartBin; ibin<=CCStopBin; ibin++) {
      binContent = max(minContent,CCSignal[ibin])*scale;
      sFit->HistEvisCCrew.at(i)->SetBinContent(ibin,binContent);
      CCInt += binContent;

      // sum up all data events
      CCdInt += sFit->HistEvisCCd.at(i)->GetBinContent(ibin);
    }

    for (Int_t ibin = NCStartBin; ibin<=NCStopBin; ibin++) {
      binContent = max( minContent,NCSignal.at(ibin) )*scale;
      sFit->HistEvisNCrew.at(i)->SetBinContent(ibin,binContent);

      if (sFit->fCutLowBins) { // cut out bins 0-2
        if (ibin>2)
          NCInt += binContent;
      }
      else
        NCInt += binContent;
      NCdInt += sFit->HistEvisNCd.at(i)->GetBinContent(ibin);
    }

    // now calculate the LL contributions (including under and overflows

    // fitting with Maximum Likelihood
    if (!sFit->fUseChi2Function) {
      MAXMSG("NCExtrapolationRS",Msg::kInfo,1) << "Using LL test..." << endl;

      // this is for fitting without the bottom two bins
      if (sFit->fCutLowBins)
        NCStartBin +=3; // cut out bins 0-2

      // contribution from NC shape
      for(Int_t ibin = NCStartBin; ibin<=NCStopBin; ibin++) {
        f -= sFit->HistEvisNCd.at(i)->GetBinContent(ibin)
          * log( sFit->HistEvisNCrew.at(i)->GetBinContent(ibin) );
      }
      // contribution from NC normalisation
      f += NCInt;
    }
    else {
      // fitting with a Chi2 function
      MAXMSG("NCExtrapolationRS",Msg::kInfo,1) << "Using Chi2 test..." << endl;
      f += sFit->GetChiSquare(0,i); // add chi2 for NC histos
    }

  }  // loop over fBeams.size()

  f *= 2; // include factor 2 in order to keep error definitions
  // as in a chi2 function
  // this also means, the UP value for errors is 1 (not 0.5)

  return;
}
//----------------------------------------------------------------------
//this method does the fit based on events in the near detector
void NCExtrapolationRS::PrepareNearDetector()
{
  //Fit happens here
  //MSG("NCExtrapolationRS",Msg::kInfo)
  //  << "Not doing fit now. Just debugging histos" << endl;

  DoFit(1);
  return;
}

//----------------------------------------------------------------------
//this method does the fit in the far detector
void NCExtrapolationRS::PredictSpectrum()
{
  PlotFinalDists();

  OscillationFit();

  return;
}
//----------------------------------------------------------------------
void NCExtrapolationRS::UseFakeData(Bool_t fake) {
  fFakeData = fake;
  if (fake) {
    for (UInt_t i=0; i<fDataFraction.size(); ++i){
      fDataFraction[BeamType::ToZarko(fBeams[i].GetBeamType())] = 1.0;
    }

    MSG("NCExtrapolationRS",Msg::kInfo) << "Running on fake data" << endl
                                        << "Setting data fraction to 1.0"
                                        << endl;
    *logFile << "Running on fake data." << endl;
  }
  else
    *logFile << "Running on REAL data." << endl;

  return;
}

//----------------------------------------------------------------------
//this method writes out fit results
void NCExtrapolationRS::WriteResources()
{
  // scale default MC histograms by bin width in order to compare to others
  for(uint i = 0; i < fBeams.size(); ++i){
    // only deal with FD
    if (fBeams[i].GetDetector() != Detector::kFar) continue;

    TH1 *defCC = fBeams[i].GetDefaultMCHistogram(NCType::kCC);
    for (int bin = 1; bin <= defCC->GetNbinsX(); ++bin) {
      Double_t val = defCC->GetBinContent(bin);
      Double_t err = defCC->GetBinError(bin);
      Double_t width = defCC->GetBinWidth(bin);
      defCC->SetBinContent(bin, val/width);
      defCC->SetBinError(bin, err/width);
    } // CC bins
    TH1 *defNC = fBeams[i].GetDefaultMCHistogram(NCType::kNC);
    for (int bin = 1; bin <= defNC->GetNbinsX(); ++bin) {
      Double_t val = defNC->GetBinContent(bin);
      Double_t err = defNC->GetBinError(bin);
      Double_t width = defNC->GetBinWidth(bin);
      defNC->SetBinContent(bin, val/width);
      defNC->SetBinError(bin, err/width);
    }  // NC bins
  } // fBeams

  NCExtrapolation::WriteResources();

  MSG("NCExtrapolationRS", Msg::kInfo)
    << "NCExtrapolationRS::WriteResources()" << endl;

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;
  TString plotsDirName = "plots"
    +NCType::kExtrapolationNames[fExtrapolationType];

  TDirectory* plotsDir = (TDirectory*) gROOT->FindObjectAny(plotsDirName.Data());
  if (plotsDir)
    plotsDir->cd();
  else
    MSG("NCExtrapolationRS", Msg::kError) << "Cannot find directory "
                                          << plotsDirName << endl;

  for(uint i = 0; i < fDataFraction.size(); ++i){
    HistEnu.at(i)->Write();
    HistEvisCC.at(i)->Write();
    HistEvisCCwt.at(i)->Write();
    HistEvisCCrew.at(i)->Write();
    HistEvisCCb.at(i)->Write();
    HistEvisCCd.at(i)->Write();
    HistEvisNC.at(i)->Write();
    HistEvisNCwt.at(i)->Write();
    HistEvisNCrew.at(i)->Write();
    HistEvisNCb.at(i)->Write();
    HistEvisNCd.at(i)->Write();

    for (uint bin = 0; bin < HistEvisCCtrueBin.at(i).size(); ++bin) {
      HistEvisCCtrueBin.at(i).at(bin)->Write();
      HistEvisNCtrueBin.at(i).at(bin)->Write();
    }
    string cName = (string) Form("cSpec_%i",i);
    TCanvas* canv =
      dynamic_cast<TCanvas*>( gROOT->FindObjectAny(cName.c_str()) );
    if (canv) canv->Write();
    else MSG("NCExtrapolationRS",Msg::kWarning) << "Can't find canvas "
                                                << cName.c_str() << endl;
    cName = (string) Form("cRat_%i",i);
    canv =
      dynamic_cast<TCanvas*>( gROOT->FindObjectAny(cName.c_str()) );
    if (canv) canv->Write();
    else MSG("NCExtrapolationRS",Msg::kWarning) << "Can't find canvas "
                                                << cName.c_str() << endl;
    cName = (string) Form("cFinalDists_%i",i);
    canv =
      dynamic_cast<TCanvas*>( gROOT->FindObjectAny(cName.c_str()) );
    if (canv) canv->Write();
    else MSG("NCExtrapolationRS",Msg::kWarning) << "Can't find canvas "
                                                << cName.c_str() << endl;


  }
  saveDir->cd();

  return;

}
//----------------------------------------------------------------------
// this method makes the coloured contributions-to-the-spectrum plots
void NCExtrapolationRS::PlotFinalDists() {
  MSG("NCExtrapolationRS",Msg::kDebug) << "Running PlotFinalDists()..."
                                       << endl;

  for(uint i = 0; i < fBeams.size(); ++i){

    // only plot stuff for ND for the moment
    if (fBeams.at(i).GetDetector()!= Detector::kNear)
      continue;

    std::vector<MCEvent>::iterator mcIt;
    for (mcIt=sFit->fMCListNear.at(i)->begin();
         mcIt<sFit->fMCListNear.at(i)->end(); mcIt++) {

      MCEvent event = (*mcIt);
      // use four wide bins in true E
      Int_t trueBin = -1;
      if (event.Enu < 4)
        trueBin = 0;
      else if (event.Enu < 8)
        trueBin = 1;
      else if (event.Enu < 15)
        trueBin = 2;
      else if (event.Enu < 120)
        trueBin = 3;
      else
        trueBin = 99;

      // include SKZP reweighting
      Double_t weight = event.BeamWeight*event.Weight;

      // do the scaling on true energy here
      if (event.McType==0 && trueBin<5)
        weight *= mScaleNC.at(trueBin);

      if ( event.SelType==0 ) {
        // this is NC selected

        HistEvisNCtrueBin.at(i).at(trueBin) -> Fill(event.Evis,weight);
      }
      else if (event.SelType==1){
        // this in CC selected

        HistEvisCCtrueBin.at(i).at(trueBin) -> Fill(event.EvisCC,weight);
      }
    } // event vector

    string canvName = (string)Form("cFinalDists_%i",i);
    TCanvas* cFinal = new TCanvas(canvName.c_str(),canvName.c_str(),
                                  40,40,500,700);
    cFinal->Divide(1,2);
    cFinal->cd(1);
    gPad->SetLogy();

    // scale according to fDataFraction and draw
    Double_t dataFrac = fDataFraction[BeamType::ToZarko(
                                        fBeams[i].GetBeamType())];
    HistEvisCCwt.at(i)->SetLineColor(2);
    if (fScaleToData.at(i)) HistEvisCCwt.at(i)->Scale(dataFrac);
    HistEvisCCwt.at(i)->Draw();
    HistEvisCCwt.at(i)->SetTitle(";E_{vis} (GeV);Events");
    HistEvisCCrew.at(i)->Draw("same");
    HistEvisCCd.at(i)->Draw("samepe");
    cFinal->cd(2);
    gPad->SetLogy();
    HistEvisNCwt.at(i)->SetLineColor(2);
    if (fScaleToData.at(i)) HistEvisNCwt.at(i)->Scale(dataFrac);
    HistEvisNCwt.at(i)->Draw();
    HistEvisNCwt.at(i)->SetTitle(";E_{vis} (GeV);Events");
    HistEvisNCrew.at(i)->Draw("same");
    HistEvisNCd.at(i)->Draw("samepe");

    TLegend *leg = new TLegend(0.6,0.6,0.9,0.9);
    leg->SetBorderSize(0);

    for (uint bin = 0; bin < HistEvisCCtrueBin.at(i).size();
         ++bin) {
      if (fScaleToData.at(i)) {
        HistEvisNCtrueBin.at(i).at(bin)->Scale(dataFrac);
        HistEvisCCtrueBin.at(i).at(bin)->Scale(dataFrac);
      }
      HistEvisNCtrueBin.at(i).at(bin)->SetLineColor(2*bin+2);
      HistEvisCCtrueBin.at(i).at(bin)->SetLineColor(2*bin+2);

      cFinal->cd(1);
      HistEvisCCtrueBin.at(i).at(bin)->Draw("same");
      cFinal->cd(2);
      HistEvisNCtrueBin.at(i).at(bin)->Draw("same");
    }

    leg->AddEntry(HistEvisNCtrueBin.at(i).at(0),"0-4  GeV","l");
    leg->AddEntry(HistEvisNCtrueBin.at(i).at(1),"4-8  GeV","l");
    leg->AddEntry(HistEvisNCtrueBin.at(i).at(2),"8-15 GeV","l");
    leg->AddEntry(HistEvisNCtrueBin.at(i).at(3),">15  GeV","l");
    leg->Draw();

  } //NBEAMS

  MSG("NCExtrapolationRS",Msg::kDebug) << "Finished PlotFinalDists()"
                                       << endl;
}

//----------------------------------------------------------------------
// this method does a grid-search oscillation fit, using the predicted
// spectrum
void NCExtrapolationRS::OscillationFit() {

  MSG("NCExtrapolationRS",Msg::kDebug) << "Running OscillationFit()..."
                                       << endl;

  // check whether we have FD beam
  Bool_t farBeamPresent = false;
  Int_t farIndex = -1;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetDetector() == Detector::kFar) {
      farBeamPresent = true;
      farIndex = i;
    }
  }

  if (!farBeamPresent) {
    MSG("NCExtrapolationRS",Msg::kError)
      << "No FD beam; not doing oscillation fit." << endl;
    return;

  }

  double chiSqr = 1.e12;

  double minChiSqr = 1.e12;
  double bestDeltaMSqr = 0.;
  double bestUMu3Sqr2Theta = 0.;
  double bestUS3Sqr = 0.;

  //get the values on the boundary
  double minBoundChiSqr = 1.e12;
  double bestBoundSterile = 0.;
  double bestBoundDeltaMSqr = 0.;

  MSG("NCExtrapolationRS", Msg::kInfo) << "NCType::kNumDeltaMSqrBins = "
                                       << NCType::kNumDeltaMSqrBins << endl;
  MSG("NCExtrapolationRS", Msg::kInfo) << "NCType::kNumUS3SqrBins = "
                                       << NCType::kNumUS3SqrBins << endl;
  MSG("NCExtrapolationRS", Msg::kInfo) << "NCType::kNumUMu3SqrBins = "
                                       << NCType::kNumUMu3SqrBins << endl;

  //store the chi^2 values in delta m^2, fs
  //to make the contour you need to find delta chi^2 for the pairs of delta m^2
  //and fs.  minimize for the systematic parameters and sin^2 2theta at each
  //pair in delta m^2 - fs space - see bevington for justification

  vector <vector <vector <double> > > chiSqrPnts;
  chiSqrPnts.resize(NCType::kNumDeltaMSqrBins);
  for(int i=0; i<NCType::kNumDeltaMSqrBins; i++) {
    chiSqrPnts[i].resize(NCType::kNumUS3SqrBins);
    for(int j=0; j < NCType::kNumUS3SqrBins; j++) {
      chiSqrPnts[i][j].resize(NCType::kNumUMu3SqrBins);
      for(int k=0; k<NCType::kNumUMu3SqrBins; k++) {
        chiSqrPnts[i][j][k] = 0.0;
      }
    }
  }

  double dmPnts[NCType::kNumDeltaMSqrBins] = {0.};
  double fsPnts[NCType::kNumUS3SqrBins] = {0.};
  double ssPnts[NCType::kNumUMu3SqrBins] = {0.};

  double DeltaMSqrWidth   = (NCType::kDeltaMSqrEnd - NCType::kDeltaMSqrStart)
    / double(NCType::kNumDeltaMSqrBins);
  double UMu3SqrWidth      = (NCType::kUMu3SqrEnd - NCType::kUMu3SqrStart)
    / double(NCType::kNumUMu3SqrBins);
  double US3SqrWidth    = (NCType::kUS3SqrEnd - NCType::kUS3SqrStart)
    / double(NCType::kNumUS3SqrBins);

  // grid-search loops
  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; ++dm) {
    dmPnts[dm] = (NCType::kDeltaMSqrStart + (float(dm+1)*DeltaMSqrWidth)
                  - DeltaMSqrWidth/2.0) * 1.0e-3;
    MSG("NCExtrapolationRS",Msg::kVerbose) << "dm = " << dmPnts[dm]
                                           << " chiSqr = " << chiSqr
                                           << " minChiSqr = " << minChiSqr
                                           << endl;

    for(int fs = 0; fs < NCType::kNumUS3SqrBins; ++fs) {
      fsPnts[fs] = NCType::kUS3SqrStart + (float(fs+1)*US3SqrWidth
                                             - US3SqrWidth/2.0);
      MSG("NCExtrapolationRS",Msg::kVerbose) << "   fs = " << fsPnts[fs]
                                             << " chiSqr = " << chiSqr
                                             << " minChiSqr = " << minChiSqr
                                             << endl;

      for(int ss = 0; ss < NCType::kNumUMu3SqrBins; ++ss){
        ssPnts[ss] = NCType::kUMu3SqrStart + (float(ss+1)*UMu3SqrWidth
                                             - UMu3SqrWidth/2.0);
        MSG("NCExtrapolationRS",Msg::kVerbose) << "      ss = " << ssPnts[ss]
                                               << " chiSqr = " << chiSqr
                                               << " minChiSqr = " << minChiSqr
                                               << endl;

        // put in oscillation parameters for energy bins here
//      fBeams[farIndex].SetOscillationParameters(ssPnts[ss],dmPnts[dm],
//                                                fsPnts[fs]);

        chiSqr = this->GetChiSqr(fBeams[farIndex]);
        chiSqrPnts[dm][fs][ss] = chiSqr;

        if(chiSqr < minChiSqr){
          minChiSqr = chiSqr;
          bestDeltaMSqr = dmPnts[dm];
          bestUMu3Sqr2Theta = ssPnts[ss];
          bestUS3Sqr = fsPnts[fs];
        }

        //get the best fit on the boundary
        if(fs == 0 && chiSqr < minBoundChiSqr){
          minBoundChiSqr = chiSqr;
          bestBoundSterile =  fsPnts[fs];
          bestBoundDeltaMSqr = dmPnts[dm];
        }

      }//end loop over sin^2 values
    }//end loop over sterile fractions
  }//end loop over delta m^2

  // set NCExtrapolation data members to best fit values
  // and oscillation parameters for NCEnergyBins to best fit
  BestDeltaMSqr() = bestDeltaMSqr;
  BestUMu3Sqr() = bestUMu3Sqr2Theta;
  BestUS3Sqr() = bestUS3Sqr;
//   fBeams[farIndex].SetOscillationParameters(fBestUMu3Sqr,fBestDeltaMSqr,
//                                          fBestUS3Sqr);

  double delta = 1.0e10;
  //double confidenceLevel = 0.;
  double minUMu3SqrChiSqr = 0., minUS3SqrChiSqr = 0.;
  int minUMu3SqrIndex = 0, minUS3SqrIndex = 0;

  //Produce Delta Chi^2 for dm2 vs fsterile
  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; ++dm){
    for(int fs = 0; fs < NCType::kNumUS3SqrBins; ++fs){
      minUMu3SqrChiSqr = 1.0e10;
      for(int ss = 0; ss < NCType::kNumUMu3SqrBins; ++ss) {
        if(chiSqrPnts[dm][fs][ss] < minUMu3SqrChiSqr) {
          minUMu3SqrChiSqr = chiSqrPnts[dm][fs][ss];
          minUMu3SqrIndex = ss;
        }
      }
      fDeltaChiSqrDeltaMVsUS3Sqr->Fill(fsPnts[fs],dmPnts[dm]*1000.0,
                                         chiSqrPnts[dm][fs][minUMu3SqrIndex]);
      if(chiSqrPnts[dm][fs][minUMu3SqrIndex] < delta) {
        delta = chiSqrPnts[dm][fs][minUMu3SqrIndex];
      }
    }
  }
  // subtract minimum chi2 -- thus getting a DeltaChi2 surface
  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; ++dm){
    for(int fs = 0; fs < NCType::kNumUS3SqrBins; ++fs){
      fDeltaChiSqrDeltaMVsUS3Sqr->SetBinContent(fs+1, dm+1,
                                              fDeltaChiSqrDeltaMVsUS3Sqr
                                                  ->GetBinContent(fs+1, dm+1)
                                                  - delta);
    }
  }

  //Produce Delta Chi^2 for dm2 vs sin2
  delta = 1.0e10;
  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; ++dm){
    for(int ss = 0; ss < NCType::kNumUMu3SqrBins; ++ss){
      minUS3SqrChiSqr = 1.0e10;
      for(int fs = 0; fs < NCType::kNumUS3SqrBins; ++fs) {
        if(chiSqrPnts[dm][fs][ss] < minUS3SqrChiSqr) {
          minUS3SqrChiSqr = chiSqrPnts[dm][fs][ss];
          minUS3SqrIndex = fs;
        }
      }
      fDeltaChiSqrDeltaMVsUMu3Sqr->Fill(ssPnts[ss], dmPnts[dm]*1000.0,
                                       chiSqrPnts[dm][minUS3SqrIndex][ss]);
      if(chiSqrPnts[dm][minUS3SqrIndex][ss] < delta) {
        delta = chiSqrPnts[dm][minUS3SqrIndex][ss];
      }
    }
  }
  // subtract minimum chi2 -- thus getting a DeltaChi2 surface
  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; ++dm){
    for(int ss = 0; ss < NCType::kNumUMu3SqrBins; ++ss){
      fDeltaChiSqrDeltaMVsUMu3Sqr->SetBinContent(ss+1, dm+1,
                                              fDeltaChiSqrDeltaMVsUMu3Sqr
                                                ->GetBinContent(ss+1, dm+1)
                                                - delta);
    }
  }


  MSG("NCExtrapolationRS", Msg::kInfo) << "best fit (f_{s},sin^2(2theta),"
                                       << "Delta m^{2}) =  ("
                                       << bestUS3Sqr
                                       << "," << bestUMu3Sqr2Theta
                                       << "," << bestDeltaMSqr
                                       << "), chi^{2} = " << minChiSqr << endl;



}

//----------------------------------------------------------------------
// this method returns the chi2 for the FD oscillation fit
double NCExtrapolationRS::GetChiSqr(NCBeam& farBeam) {

  MSG("NCExtrapolationRS",Msg::kVerbose) << "Running GetChiSqr()"
                                         << endl;

  double chiSqr = 0.;
  //  double bg = 0.;
  double mc = 0.;
  double signal = 0.;
  double signalError = 0.;

  //loop over the bins to find the chisqr
  for(int j = 0; j < farBeam.GetNumberEnergyBins(NCType::kCC); ++j) {

    NCEnergyBin* energyBin = farBeam.GetEnergyBin(j, NCType::kCC);
    signal = energyBin->GetSignal();
    signalError = energyBin->GetSignalError();

    //mc already contains the value of the background as the
    //GetMCExpectationForBin returns the sum of MC signal and
    //background
    mc = energyBin->GetTotalMC();
//     bg = energyBin->GetMCBackground() + energyBin->GetMCTau() + energyBin->GetMCElectron();

    //the following is the chi^2 for poisson distributed data:
    if(mc <= 0.) mc = 1.e-3;

    if(mc > 0.){
      chiSqr += mc;
      if(signal > 0.) chiSqr += signal*(TMath::Log( signal/mc ) - 1.);
      if (signal == 0.) chiSqr += -1*signal;
    }
  } // CC loop

  // NC loop
  for(int j = 0; j < farBeam.GetNumberEnergyBins(NCType::kNC); ++j) {

    NCEnergyBin* energyBin = farBeam.GetEnergyBin(j, NCType::kNC);
    signal = energyBin->GetSignal();
    signalError = energyBin->GetSignalError();

    //mc already contains the value of the background as the
    //GetMCExpectationForBin returns the sum of MC signal and
    //background
    mc = energyBin->GetTotalMC();
//     bg = energyBin->GetMCBackground() + energyBin->GetMCTau() + energyBin->GetMCElectron();

    //the following is the chi^2 for poisson distributed data:
    if(mc < 0.) mc = 1.e-3;

    if(mc > 0.){
      chiSqr += mc;
      if(signal > 0.) chiSqr += signal*(TMath::Log( signal/mc ) - 1.);

    }
  } // NC loop

  //have to multiply by 2
  chiSqr *= 2.;

  MSG("NCExtrapolationRS", Msg::kVerbose) << "chi^2 = " << chiSqr << endl;





/*

  //add in the penalty terms

  for(uint i = 0; i < fParameters.size(); ++i){
    if(fFitParameters[i]){
      chiSqr += (fParameters[i]*fParameters[i])/(fParameterSigmas[i]*fParameterSigmas[i]);

      MSG("NCExtrapolationRS", Msg::kDebug) << i << " "
                                            <<  fParameters[i] << " "
                                            << fParameterSigmas[i] << " "
                                            << chiSqr << endl;
    }
  }

  if(fFittingShower) chiSqr += ((kChangeValuesShower[fShowerSystematic]*kChangeValuesShower[fShowerSystematic]) /
                                (fParameterSigmas[NCType::kShowerEnergy]*fParameterSigmas[NCType::kShowerEnergy]));
  if(fFittingTrack) chiSqr += ((kChangeValuesTrack[fTrackSystematic]*kChangeValuesTrack[fTrackSystematic]) /
                                (fParameterSigmas[NCType::kTrackEnergy]*fParameterSigmas[NCType::kTrackEnergy]));
*/

  return chiSqr;
}

---