NCExtrapolationDP.cxx

Go to the documentation of this file.

//$Id: NCExtrapolationDP.cxx,v 1.35 2008/01/11 15:53:31 bckhouse Exp $
//
//NCExtrapolationDP.cxx
//
//class to hold pdfs for doing nccc separation
//
//B. Rebel 2/2007

#include "NCUtils/Extrapolation/NCExtrapolationDP.h"
#include "MessageService/MsgService.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 "TStopwatch.h"

#include <cassert>


//=================David: please add code to return the histograms/graphs
//with the predicted spectrum, background spectrum, fit spectrum, and data
//the interfaces are defined in NCExtrapolation.h
//
//the flow of this code is seen in lines NCExtrapolationModule::PrepareForExtrapolation
// and NCExtrapolationModule::ExtrapolateNearToFar
//
//the main difference here is that i try to fill as many of the histograms/arrays you
//use in the different macros at once.  so the near det arrays are filled at the same
//time and the far det ones are too.  there is still a separation between the fitting
//of the ND data and the filling of the far detector MC so that it can be adjusted for
//the returned fit values.


static const int kNorm = 0;
static const int kEshower = 1;
static const int kEmu = 2;
static const int kMA_QE = 3;
static const int kMA_RES = 4;
static const int kdisfact = 5;
static const int kBeam1 = 6;
static const int kBeam2 = 7;
static const int kBeam3 = 8;
static const int kBeam4 = 9;
static const int kBeam5 = 10;
static const int kReco = 11;
static const int knc = 12;
static TString kParNames[13] = {"Normalization", "E_shower",
                                "E_muon", "ma_qe", "ma_res",
                                "DISFACT", "Beam 1", "Beam 2",
                                "Beam 3", "Beam 4", "Beam 5",
                                "Reco", "NC"};

using namespace NCType;

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

//......................................................................
NCExtrapolationDP::NCExtrapolationDP() :
  NCExtrapolation::NCExtrapolation()
{
}

//.......................................................................
NCExtrapolationDP::NCExtrapolationDP(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)
{
  fExtrapolationType = NCType::kExtrapolationDP;
  fsys = 0;
}

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

//----------------------------------------------------------------------
//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 NCExtrapolationDP::Prepare(const Registry& r)
{

  NCExtrapolation::Prepare(r);

  int tmpi;
  double tmpd;
  if(r.Get("Systematics", tmpi)) fsys         = tmpi;
  if(r.Get("FarPOT",      tmpd)) fpot         = tmpd;
  if(r.Get("MaxFarEReco", tmpd)) fMaxFarEReco = tmpd;

  double y0[5]={0,0.25,0.5,0.75,1};

  fnom = new TH2F("maqenom","", NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear,4,y0);
  fmaqewgt = new TH2F("maqewgt","", NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear,4,y0);
  fmareswgt = new TH2F("mareswgt","", NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear,4,y0);
  fdisfactwgt = new TH2F("disfactwgt","", NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear,4,y0);

  fnomcc=new TH1F("maqenomcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmaqewgtcc=new TH1F("maqewgtcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmareswgtcc=new TH1F("mareswgtcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fdisfactwgtcc=new TH1F("disfactwgtcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);

  fnomnc=new TH1F("maqenomnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmaqewgtnc=new TH1F("maqewgtnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmareswgtnc=new TH1F("mareswgtnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fdisfactwgtnc=new TH1F("disfactwgtnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);

  fdata=new TH1D("data","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fdatanc=new TH1D("datanc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);

  fmcnomcc=new TH1D("mcnomcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmcbfcc=new TH1D("mcbfcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmtempcc=new TH1D("mtempcc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);

  fmcnomnc=new TH1D("mcnomnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmcbfnc=new TH1D("mcbfnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);
  fmtempnc=new TH1D("mtempnc","",NCType::kNumEnergyBinsNear,NCType::kEnergyBinsNear);

  for (int x=0;x<21*21;x++) {
    for (int y=0;y<NCType::kNumEnergyBinsNear;y++) {

      fncvcc[x][y]=0;
      fncvnc[x][y]=0;

      for (int z=0;z<4;z++) {

        fvexp[x][y][z] = 0.;
        fvnc[x][y][z] = 0.;

      }
    }
  }

  for (int i=0;i<NCType::kNumEnergyBinsFar;i++) {
    fobsall[0][i] = 0.;
    fobsall[1][i] = 0.;
  }

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationDP::AddEvent(ANtpHeaderInfo *headerInfo,
                                 ANtpBeamInfo *beamInfo,
                                 ANtpRecoInfo *recoInfo,
                                 ANtpAnalysisInfo *analysisInfo,
                                 ANtpTruthInfoBeam *truthInfo,
                                 int runType,
                                 bool useMCAsData,
                                 int fileType)
{
  //++++Energies are set to the correct NC/CC values in the NCExtrapolationModule
  //before being passed into the extrapolations

  //only use the L010z185i beam - taus not used yet either
  //dont bother with events outside the fiducial volume
  //or those that are too high of energy
  if(beamInfo->beamType != BeamType::ToZarko(BeamType::kL010z185i)
     || recoInfo->inFiducialVolume < 1
     || recoInfo->nuEnergy > NCType::kMaxEnergy) return;

  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 NCExtrapolationDP::AddNearEvent(ANtpHeaderInfo *headerInfo,
                                     ANtpBeamInfo *beamInfo,
                                     ANtpRecoInfo *recoInfo,
                                     ANtpAnalysisInfo *analysisInfo,
                                     ANtpTruthInfoBeam *truthInfo,
                                     int  runType,
                                     bool useMCAsData,
                                     int fileType)
{
  MSG("NCExtrapolationDP", Msg::kDebug) << fileType << " " << beamInfo->beamType << endl;

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

  BeamType::BeamType_t beamType = BeamType::FromZarko(beamInfo->beamType);
  int beam = NCExtrapolation::FindBeamIndex(Detector::kNear, beamType, runType);
  double scaleFactor = fDataPOTNear[beamInfo->beamType]/fMCPOTNear[beamInfo->beamType];

  if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
    //compute the weight for a 15% change each in ma_qe and ma_res
    //and a 20% change in the DISFACT parameters
    vector<bool> useParameters;
    vector<double> adjustedValues;

    for(int i = 0; i < NCType::kNumNeugenParameters; ++i){
      useParameters.push_back(false);
      adjustedValues.push_back(NCType::kParameterDefaults[i]);
    }

    int pars[] = {NCType::kma_qe, NCType::kma_res, NCType::kDISFACT};
    double weights[] = {1., 1., 1.};
    double shifts[] = {1.15, 1.15, 1.20};

    for(int i = 0; i < 3; ++i){
      adjustedValues[pars[i]] = shifts[i]*NCType::kParameterDefaults[pars[i]];
      useParameters[pars[i]] = true;
      weights[i] = fUtils->FindNeugenWeight(useParameters, adjustedValues);
      adjustedValues[pars[i]] = NCType::kParameterDefaults[pars[i]];
      useParameters[pars[i]] = false;
    }

    //==============David: i added the following as i cant see that your
    //code accounts for it
    //check if this is nc-like or cc-like to get the right event energy
    float ehad=0.;
    float emu=0.;

    float wgtpt=0.;
    float wgt2 = 1.;
    float prob = 1.;

    //===============David:  please make the following
    //settable by using the registry passed into the Prepare method
    //sys should also by settable using the registry
    int esmin=10;
    int esmax=11;
    int mumin=10;
    int mumax=11;
    if(fsys==1){
      esmin=0;
      esmax=21;
      mumin=0;
      mumax=21;
    }
    float ehadshift=0.;
    float emushift=0.;
    float ereco = 0.;
    float recoy = 0.;
    int b = 0;
    int ybin =0;

    if(analysisInfo->isNC > 0){

      ehad = recoInfo->showerEnergy;
      emu = recoInfo->muEnergy;
      wgtpt = recoInfo->weight;

      wgt2=1; //maqewgt*mareswgt;
      prob=1.0;
      prob*=wgtpt;

      fmcnomnc->Fill(ehad+emu,prob);

      MAXMSG("NCExtrapolationDP", Msg::kInfo, 10) << "mc " << recoInfo->nuEnergy
                                                  << "/" << recoInfo->showerEnergyNC
                                                  << "/" << recoInfo->showerEnergy
                                                  << " " << recoInfo->muEnergy
                                                  << " " << recoInfo->weight
                                                  << "/" << prob << endl;

      if(truthInfo->interactionType == NCType::kNC){
        fmaqewgtnc->Fill(ehad+emu, weights[0]);
        fmareswgtnc->Fill(ehad+emu, weights[1]);
        fdisfactwgtnc->Fill(ehad+emu, weights[2]);
        fnomnc->Fill(ehad+emu);
      }
      else if(truthInfo->interactionType == NCType::kCC){
        fmaqewgtcc->Fill(ehad+emu, weights[0]);
        fmareswgtcc->Fill(ehad+emu, weights[1]);
        fdisfactwgtcc->Fill(ehad+emu, weights[2]);
        fnomcc->Fill(ehad+emu);
      }

      ehadshift = 0.;
      emushift = 0.;

      for(int q = esmin; q < esmax; q++){

        ehadshift = ehad*(0.8 + q*0.02);

        for(int l = mumin; l < mumax; l++){

          emushift = emu*(0.9 + l*0.01);

          ereco = ehadshift+emushift;

          b = FindRecoEnergyBin(ereco, beam);

          if(truthInfo->interactionType == NCType::kCC) {
            fncvcc[q*21+l][b] += prob;
          }

          if (truthInfo->interactionType == NCType::kNC) {
            fncvnc[q*21+l][b] += prob;
          }

        }//end loop over muon energies
      }//end loop over shower energies
    }//end if nc-like
    else if(analysisInfo->isCC > 0){

      ehad = recoInfo->showerEnergy;
      if(ehad < 0.) ehad = 0.;
      emu = recoInfo->muEnergy;
      wgtpt = recoInfo->weight;

      wgt2 = 1.; //maqewgt*mareswgt;
      prob = 1.0;
      prob *= wgtpt;

      fmcnomcc->Fill(recoInfo->nuEnergy, prob);
      fmaqewgt->Fill(recoInfo->nuEnergy, recoInfo->hadronicY, weights[0]);
      fmareswgt->Fill(recoInfo->nuEnergy, recoInfo->hadronicY, weights[1]);
      fdisfactwgt->Fill(recoInfo->nuEnergy, recoInfo->hadronicY, weights[2]);
      fnom->Fill(recoInfo->nuEnergy, recoInfo->hadronicY);

      ehadshift = 0.;
      emushift = 0.;

      for(int q = esmin; q < esmax; q++){

        ehadshift = ehad*(0.8 + q*0.02);

        for(int l = mumin; l < mumax; l++){

          emushift = emu*(0.9+l*0.01);

          ereco = ehadshift + emushift;

          recoy = ehadshift/ereco;
          ybin = TMath::Nint(4*recoy);
          if(ybin > 3) ybin = 3;

          b = FindRecoEnergyBin(ereco, beam);

          fvexp[q*21+l][b][ybin] += prob*scaleFactor;

          if (truthInfo->interactionType == NCType::kNC) {
            fvnc[q*21+l][b][ybin] += prob*scaleFactor;
          }
        }//end loop over muon energies
      }//end loop over shower energies
    }//end if cc-like
  }//end if MC
  else if(headerInfo->dataType == (int)SimFlag::kData || useMCAsData){

    //when using MC as data the module sets the recoInfo->weight to be
    //the weight for the input value.  when using data recoInfo->weight = 1.
    //FinalEventCheck in module sets recoInfo->showerEnergy, muEnergy and nuEnergy
    //to correct values for nc or cc
    if(analysisInfo->isNC > 0){
      fdatanc->Fill(recoInfo->nuEnergy, recoInfo->weight);
      MAXMSG("NCExtrapolationDP", Msg::kInfo, 10) << "data " << recoInfo->nuEnergy
                                                  << "/" << recoInfo->showerEnergyNC
                                                  << "/" << recoInfo->showerEnergy
                                                  << " " << recoInfo->muEnergy
                                                  << " " << recoInfo->weight
                                                  << endl;
    }
    else if(analysisInfo->isCC > 0){
      fdata->Fill(recoInfo->nuEnergy, recoInfo->weight);
    }

  }//end if data

  return;
}

//----------------------------------------------------------------------
void NCExtrapolationDP::AddFarEvent(ANtpHeaderInfo *headerInfo,
                                    ANtpBeamInfo *beamInfo,
                                    ANtpRecoInfo *recoInfo,
                                    ANtpAnalysisInfo *analysisInfo,
                                    ANtpTruthInfoBeam *truthInfo,
                                    int runType,
                                    bool useMCAsData,
                                    int fileType)
{
  MAXMSG("NCExtrapolationDP", Msg::kInfo, 1) << "AddFarEvent " << endl;
  MSG("NCExtrapolationDP", Msg::kDebug) << fileType << " " << beamInfo->beamType << endl;

  //need to check if this type of beam already exists in the vector.
  //if so add the event to that beam. if not, make it the beam then
  //the add event to it

  //check on data/mc, if mc adjust the weights/energies to reflect fit in nd
  if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){

    //reweight the event according to the ND fit
    vector<bool> useParameters;
    vector<double> adjustedValues;

    int pars[] = {NCType::kma_qe, NCType::kma_res, NCType::kDISFACT};
    double shifts[] = {fNDFitResults[kMA_QE]*0.2 - 4.,
                       fNDFitResults[kMA_RES]*0.2 - 4.,
                       fNDFitResults[kdisfact]*0.2 - 4.};

    for(int i = 0; i < NCType::kNumNeugenParameters; ++i){
      useParameters.push_back(false);
      adjustedValues.push_back(NCType::kParameterDefaults[i]);
    }

    for(int i = 0; i < 3; ++i){
      useParameters[pars[i]] = true;
      adjustedValues[pars[i]] = shifts[i]*NCType::kParameterSigmas[pars[i]];
      adjustedValues[pars[i]] += NCType::kParameterDefaults[pars[i]];
      MAXMSG("NCExtrapolationDP", Msg::kInfo, 3) << adjustedValues[pars[i]] << " "
                                                 << useParameters[pars[i]] << " "
                                                 << NCType::kParameterNames[pars[i]]
                                                 << endl;
    }

    double weightNeugen = fUtils->FindNeugenWeight(useParameters, adjustedValues);
    double weightNC = 1.;
    if(truthInfo->interactionType == NCType::kNC)
      weightNC = 0.8+0.02*fNDFitResults[knc];
    double weightNorm = 0.9 + 0.005*fNDFitResults[kNorm];

    //recoInfo->weight has the SKZP weight for the event
    double totalWeight = weightNeugen*weightNC*weightNorm*recoInfo->weight;

    //FinalEventCheck in module sets recoInfo->showerEnergy, muEnergy and nuEnergy
    //to correct values for nc or cc
    //double nominalEnergy = recoInfo->nuEnergy;
    recoInfo->showerEnergy *= 0.8 + 0.02*fNDFitResults[kEshower];
    recoInfo->muEnergy *= 0.9 + 0.01*fNDFitResults[kEmu];
    recoInfo->nuEnergy = recoInfo->muEnergy + recoInfo->showerEnergy;

    MAXMSG("NCExtrapolationDP", Msg::kInfo,1) << 0.8 + 0.02*fNDFitResults[kEshower] << " "
                                              << 0.9 + 0.01*fNDFitResults[kEmu] << " "
                                              << 0.9 + 0.005*fNDFitResults[kNorm] << endl;
    MAXMSG("NCExtrapolationDP", Msg::kInfo,1) << " weight NC = " << weightNC
                                              << " weight norm = " << weightNorm
                                              << " weight neugen = " << weightNeugen
                                              << " recoInfo->weight = " << recoInfo->weight
                                              << " total weight = " << totalWeight
                                              << endl;

  }//end if mc

  if(recoInfo->nuEnergy > fMaxFarEReco){
    MAXMSG("NCExtrapolationDP", Msg::kInfo, 100) << "****Max far E_reco exceeded :"
                                                 << recoInfo->nuEnergy << "/"
                                                 << fMaxFarEReco << endl;
    return;
  }

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

  return;
}

//----------------------------------------------------------------------
//this method does the fit based on events in the near detector
void NCExtrapolationDP::PrepareNearDetector()
{

  //find the position in the beam vector, then make the data graph
  int runType = 1;
  int beam = NCExtrapolation::FindBeamIndex(Detector::kNear, BeamType::kL010z185i, runType);
  if(beam > 100){
    MSG("NCExtrapolationDP", Msg::kWarning) << "could not find "
                                            << "ND beam in vector - bail"
                                            << endl;
    return;
  }

  //divide the histograms for the weights
  fmaqewgt->Divide(fnom);
  fmareswgt->Divide(fnom);
  fdisfactwgt->Divide(fnom);

  fmaqewgtcc->Divide(fnomcc);
  fmareswgtcc->Divide(fnomcc);
  fdisfactwgtcc->Divide(fnomcc);

  fmaqewgtnc->Divide(fnomnc);
  fmareswgtnc->Divide(fnomnc);
  fdisfactwgtnc->Divide(fnomnc);

  //=========David:  this code comes from ndfit.C. it would be
  //good if you could add some comments about what you are doing here.
  float qewgt[NCType::kNumEnergyBinsNear][4];
  float reswgt[NCType::kNumEnergyBinsNear][4];
  float diswgt[NCType::kNumEnergyBinsNear][4];


  for (int i=0;i<NCType::kNumEnergyBinsNear;i++) {
    for (int j=0;j<4;j++) {
      qewgt[i][j]=fmaqewgt->GetBinContent(i+1,j+1)-1.0;
      reswgt[i][j]=fmareswgt->GetBinContent(i+1,j+1)-1.0;
      diswgt[i][j]=fdisfactwgt->GetBinContent(i+1,j+1)-1.0;
    }
  }


  float qewgtnc_cc[NCType::kNumEnergyBinsNear];
  float reswgtnc_cc[NCType::kNumEnergyBinsNear];
  float diswgtnc_cc[NCType::kNumEnergyBinsNear];

  float qewgtnc_nc[NCType::kNumEnergyBinsNear];
  float reswgtnc_nc[NCType::kNumEnergyBinsNear];
  float diswgtnc_nc[NCType::kNumEnergyBinsNear];

  for(int i = 0; i < NCType::kNumEnergyBinsNear; i++){

    qewgtnc_cc[i] = fmaqewgtcc->GetBinContent(i+1)-1.0;
    reswgtnc_cc[i] = fmareswgtcc->GetBinContent(i+1)-1.0;
    diswgtnc_cc[i] = fdisfactwgtcc->GetBinContent(i+1)-1.0;

    qewgtnc_nc[i] = fmaqewgtnc->GetBinContent(i+1)-1.0;
    reswgtnc_nc[i] = fmareswgtnc->GetBinContent(i+1)-1.0;
    diswgtnc_nc[i] = fdisfactwgtnc->GetBinContent(i+1)-1.0;

  }

  float ccefficarr[11][NCType::kNumEnergyBinsNear];
  float ncefficarr[11][NCType::kNumEnergyBinsNear];
  float ccpurarr[11][NCType::kNumEnergyBinsNear];

  for(int i = 0; i < 11; i++){
    for(int j = 0; j < NCType::kNumEnergyBinsNear; j++){
      ccefficarr[i][j] = 1.0;
      ncefficarr[i][j] = 1.0;
      ccpurarr[i][j] = 1.0;
    }
  }

  // data/MC normalisation factor
  float factle=1.0; // le

  // evis distribution for "data"
  // le-10
  for(int j = 0; j < NCType::kNumEnergyBinsNear; j++){
    fvobs[j]=fdata->GetBinContent(j+1);
    fvobsnclike[j]=fdatanc->GetBinContent(j+1);
  }

  // define evis histograms for "MC" sample
  //le-10 - cc

  float chisqmine=9999;
  float chisqmincc=9999;
  float chisqminnc=9999;
  float chisq=0;
  float chisqnc=0;

  int cnorm=0;
  int cesh=0;
  int cemu=0;
  int cmaqe=0;
  int cmares=0;
  int cdis=0;
  int cbeam1=0;
  int cbeam2=0;
  int cbeam3=0;
  int cbeam4=0;
  int cbeam5=0;

  int creco=0;
  int cnc=0;

  float beamerr[NCType::kNumEnergyBinsNear][5];

  for(int i = 0; i < NCType::kNumEnergyBinsNear; i++){
    beamerr[i][0] = 1.0;
    beamerr[i][1] = 1.0;
    beamerr[i][2] = 1.0;
    beamerr[i][3] = 1.0;
    beamerr[i][4] = 1.0;
  }

  int nomin = 20;
  int nomax = 21;
  int esmin = 10;
  int esmax = 11;
  int mumin = 10;
  int mumax = 11;
  int dmmin = 20;
  int dmmax = 21;
  int dismin = 20;
  int dismax = 21;

  if(fsys==1){
    nomin = 0;
    nomax = 41;
    esmin = 0;
    esmax = 21;
    mumin = 0;
    mumax = 21;
    dmmin = 0;
    dmmax = 41;
    dismin = 20;
    dismax = 41;
  }

  float normfact = 0.;
  float eshift = 0.;
  float emushift = 0.;
  float maqe = 0.;
  float mares = 0.;
  float disfact = 0.;
  float beamunc1 = 0.;
  float beamunc2 = 0.;
  float beamunc3 = 0.;
  float beamunc4 = 0.;
  float beamunc5 = 0.;
  float ncunc = 0.;
  int recerr = 5; //makes life easier

  float so = 0.;
  float se = 0.;
  float snc = 0.;
  float obs = 0.;
  float ex = 0.;
  float expnc = 0.;
  float beamfact = 0.;
  float ncex = 0.;
  float ncexpnc = 0.;
  float sonc = 0.;
  float senc = 0.;
  float chisqall = 0.;
  float cceff = 0.;
  float ccpur = 0.;
  float nceff = 0.;

  // loop 1 - OVERALL NORMALISATION
  for(int no = nomin; no < nomax; no++){
    normfact = 0.9 + 0.005*no;
    MAXMSG("NCExtrapolationDP", Msg::kDebug, nomax-nomin)
      << "normfact " << normfact << "/" << no << endl;

    // loop 4 - MAQEL
    for(int dm = dmmin; dm < dmmax; dm++){
      maqe = -4. + 0.2*dm;
      mares = -4. + 0.2*dm;
      MAXMSG("NCExtrapolationDP", Msg::kDebug, dmmax-dmmin)
        << "maqe,mares " << maqe << "/" << dm << endl;

      // loop 5 - DISFACT
      for(int dis = dismin; dis < dismax; dis++){
        disfact = -4 + 0.2*dis;
        MAXMSG("NCExtrapolationDP", Msg::kDebug, dismax-dismin)
          << "disfact " << disfact << "/" << dis << endl;

        // loop 2 - HADRON ENERGY SCALE
        for(int es = esmin; es < esmax; es++){
          eshift = es*0.02 + 0.8;

          // loop 3 - MUON ENERGY SCALE
          for(int em = mumin; em < mumax; em++){
            emushift = em*0.01 + 0.9;

            fmtempcc->Reset();
            fmtempnc->Reset();

            // CALCULATE CHISQ
            chisq = 0.;
            chisqnc = 0.;

            // CC-LIKE
            for(int i = 0; i < NCType::kNumEnergyBinsNear; i++){

              // le beamfact
              beamfact = (1+beamunc1*(beamerr[i][0]-1));
              beamfact *= (1+beamunc2*(beamerr[i][1]-1));
              beamfact *= (1+beamunc3*(beamerr[i][2]-1));
              beamfact *= (1+beamunc4*(beamerr[i][3]-1));
              beamfact *= (1+beamunc5*(beamerr[i][4]-1));

              // le counters
              so = 0.;
              se = 0.;
              snc = 0.;
              sonc = 0.;
              senc = 0.;

              // cc-like
              obs = fvobs[i];
              so += obs;

              // nc-like
              sonc = fvobsnclike[i];

              cceff = ccefficarr[recerr][i];
              ccpur = ccpurarr[recerr][i];
              nceff = ncefficarr[recerr][i];

              // Y-loop
              for(int j = 0; j < 4; j++){
                // calculate expectation - LE-10

                ex = fvexp[es*21+em][i][j]*normfact*beamfact*factle;
                ex *= (1.0 + qewgt[i][j]*maqe);
                ex *= (1.0 + reswgt[i][j]*mares);
                ex *= (1.0 + diswgt[i][j]*disfact);

                // apply effic/pur correction
                ex *= cceff/ccpur;

                expnc = fvnc[es*21+em][i][j]*normfact*beamfact*factle;

                // apply nc purity
                expnc *= nceff;

                //expeced cc signal
                se += ex;

                //expected nc bg to cc
                snc += expnc*ncunc;

              }//end loop over y bins

              // nc
              //expected cc bg to nc signal
              ncex = fncvcc[es*21+em][i]*normfact*beamfact*factle;
              ncex *= (1.0 + qewgtnc_cc[i]*maqe);
              ncex *= (1.0 + reswgtnc_cc[i]*mares);
              ncex *= (1.0 + diswgtnc_cc[i]*disfact);

              //expected nc signal
              ncexpnc=fncvnc[es*21+em][i]*normfact*beamfact*factle;
              ncexpnc *= (1.0 + qewgtnc_nc[i]*maqe);
              ncexpnc *= (1.0 + reswgtnc_nc[i]*mares);
              ncexpnc *= (1.0 + diswgtnc_nc[i]*disfact);
              ncexpnc *= 1.+ncunc;
              senc = ncex + ncexpnc;

              fmtempcc->SetBinContent(i+1,se+snc);
              fmtempnc->SetBinContent(i+1,senc);

              // le chisq
              if(se + snc > 0){
                chisq += TMath::Power(so-(se+snc),2)/(se+snc);
              }

              if(senc > 0){
                chisqnc += TMath::Power(sonc-senc,2)/(senc);
              }

            }//end loop over energy bins

            chisq += TMath::Power(maqe,2);
            chisq += TMath::Power(mares,2);
            chisq += TMath::Power(disfact,2);
            chisq += TMath::Power(eshift-1,2)/TMath::Power(0.10,2);
            chisq += TMath::Power(emushift-1,2)/TMath::Power(0.02,2);
            chisq += TMath::Power(normfact-1,2)/TMath::Power(0.05,2);

            chisqall = chisq + chisqnc;

            if(chisqall < chisqmine){

              chisqmine = chisqall;
              chisqmincc = chisq;
              chisqminnc = chisqnc;

              cnorm = no;
              cesh = es;
              cemu = em;
              cmaqe = dm;
              cmares = dm;
              cdis = dis;
              cbeam1 = 10;
              cbeam2 = 10;
              cbeam3 = 10;
              cbeam4 = 10;
              cbeam5 = 10;

              creco = 5;
              cnc = 10;
              fmcbfcc->Reset();
              fmcbfnc->Reset();

              for(int ll = 0; ll < NCType::kNumEnergyBinsNear; ll++){
                fmcbfcc->SetBinContent(ll+1, fmtempcc->GetBinContent(ll+1));
                fmcbfnc->SetBinContent(ll+1, fmtempnc->GetBinContent(ll+1));
              }//end loop over energy bins to fill output histogram
            }//end if min chi^2
          } // end of dis loop
        } // end of maqe + mares loop
      } // end of emu loop
    } // end of ehad loop
  } // end of norm loop

  fNDFitResults.clear();
  fNDFitResults.push_back(cnorm);
  fNDFitResults.push_back(cesh);
  fNDFitResults.push_back(cemu);
  fNDFitResults.push_back(cmaqe);
  fNDFitResults.push_back(cmares);
  fNDFitResults.push_back(cdis);
  fNDFitResults.push_back(cbeam1);
  fNDFitResults.push_back(cbeam2);
  fNDFitResults.push_back(cbeam3);
  fNDFitResults.push_back(cbeam4);
  fNDFitResults.push_back(cbeam5);
  fNDFitResults.push_back(creco);
  fNDFitResults.push_back(cnc);

  MSG("NCExtrapolationDP", Msg::kInfo)
    << "--------Results of NDFit---------"
    << endl << "chisqrmin = " << chisqmine
    << " cc-chisqr = " << chisqmincc
    << " nc-chisqr = " << chisqminnc
    << endl << kParNames[kNorm] << " = " << fNDFitResults[kNorm]*0.005+0.9 - 1.
    << endl << kParNames[kEshower] << " = "
    << (fNDFitResults[kEshower]*0.02+0.8-1.)*NCType::kParameterSigmas[NCType::kShowerEnergy]
    +NCType::kParameterDefaults[NCType::kShowerEnergy]
    << " = " << (fNDFitResults[kEshower]*0.02+0.8)
    << endl << kParNames[kEmu] << " = "
    << (fNDFitResults[kEmu]*0.01+0.9-1.)*NCType::kParameterSigmas[NCType::kTrackEnergy]
    +NCType::kParameterDefaults[NCType::kTrackEnergy]
    << " = " << (fNDFitResults[kEmu]*0.01+0.9)
    << endl << kParNames[kMA_QE] << " = "
    << (fNDFitResults[kMA_QE]*0.2-4.)*NCType::kParameterSigmas[NCType::kma_qe]
    + NCType::kParameterDefaults[NCType::kma_qe]
    << " = " << (fNDFitResults[kMA_QE]*0.2-4.)
    << endl << kParNames[kMA_RES] << " = "
    << (fNDFitResults[kMA_RES]*0.2-4.)*NCType::kParameterSigmas[NCType::kma_res]
    + NCType::kParameterDefaults[NCType::kma_res]
    << " = " << (fNDFitResults[kMA_RES]*0.2-4.)
    << endl << kParNames[kdisfact] << " = "
    << (fNDFitResults[kdisfact]*0.2-4.)*NCType::kParameterSigmas[NCType::kDISFACT]
    + NCType::kParameterDefaults[NCType::kDISFACT]
    << " = " << (fNDFitResults[kdisfact]*0.2-4.)
    << endl << kParNames[kBeam1] << " = "
    << (fNDFitResults[kBeam1]*0.4-4.)*1.
    << endl << kParNames[kBeam2] << " = "
    << (fNDFitResults[kBeam2]*0.4-4.)*1.
    << endl << kParNames[kBeam3] << " = "
    << (fNDFitResults[kBeam3]*0.4-4.)*1.
    << endl << kParNames[kBeam4] << " = "
    << (fNDFitResults[kBeam4]*0.4-4.)*1.
    << endl << kParNames[kBeam5] << " = "
    << (fNDFitResults[kBeam5]*0.4-4.)*1.
    << endl << kParNames[kReco] << " = "
    << (fNDFitResults[kReco]*0.4-2.)*1.
    << endl << kParNames[knc] << " = "
    << (fNDFitResults[knc]*0.02-0.2)*NCType::kParameterSigmas[NCType::kNCBackground]
    + NCType::kParameterDefaults[NCType::kNCBackground]
    << endl;

  //loop over the bins one more time to set the fit values in the beams
  NCEnergyBin *bin = 0;
  for(int i = 0; i < NCType::kNumEnergyBinsNear; ++i){

    // le counters
    se = 0;
    snc = 0;
    ex = 0;
    expnc = 0;
    ncex = 0;
    ncexpnc = 0;

    cceff = ccefficarr[creco][i];
    ccpur = ccpurarr[creco][i];
    nceff = ncefficarr[creco][i];

    // Y-loop
    for(int j = 0; j < 4; j++){
      // calculate expectation - LE-10
      ex = fvexp[cesh*21+cemu][i][j]*(cnorm*0.005+0.9)*factle;
      ex *= (1.0+qewgt[i][j]*(-4.+cmaqe*0.2));
      ex *= (1.0+reswgt[i][j]*(-4.+cmares*0.2));
      ex *= (1.0+diswgt[i][j]*(-4.+cdis*0.2));

      // apply effic/pur correction
      ex *= cceff/ccpur;

      //expected nc background
      expnc = fvnc[cesh*21+cemu][i][j]*(cnorm*0.005+0.9)*factle;

      // apply nc purity
      expnc *= nceff;

      //expeced cc signal
      se += ex;

      //expected nc bg to cc
      snc += expnc*(-0.2+0.02*cnc);

    }

    // nc
    //expected cc bg to nc signal
    ncex = fncvcc[cesh*21+cemu][i]*(cnorm*0.005+0.9)*factle;
    ncex *= (1.0+qewgtnc_cc[i]*(-4.+0.2*cmaqe));
    ncex *= (1.0+reswgtnc_cc[i]*(-4.+0.2*cmares));
    ncex *= (1.0+diswgtnc_cc[i]*(-4.+0.2*cdis));

    //expected nc signal
    ncexpnc = fncvnc[cesh*21+cemu][i]*(cnorm*0.005+0.9)*factle;
    ncexpnc *= (1.0+qewgtnc_nc[i]*(-4.+0.2*cmaqe));
    ncexpnc *= (1.0+reswgtnc_nc[i]*(-4.+0.2*cmares));
    ncexpnc *= (1.0+diswgtnc_nc[i]*(-4.+0.2*cdis));
    ncexpnc *= 1.+(-0.2+0.02*cnc);

    //electrons are accounted for in totals already
    bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
//     bin->SetMCFit(se);
//     bin->SetMCFitBackground(snc);
//     bin->SetMCFitElectron(0.);

    bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
//     bin->SetMCFit(ncexpnc);
//     bin->SetMCFitBackground(ncex);
//     bin->SetMCFitElectron(0.);


  }//end loop over energy bins to set fit and fit background

  return;
}

//----------------------------------------------------------------------
//this method does the fit in the far detector
void NCExtrapolationDP::PredictSpectrum()
{
  MSG("NCExtrapolationDP", Msg::kWarning) << "Predict Spectrum"
                                          << endl;

  //find the position in the beam vector
  int runType = 1;
  int beam = NCExtrapolation::FindBeamIndex(Detector::kFar, BeamType::kL010z185i, runType);

  //uncomment next 3 lines if you dont want to fit for oscillation parameters
  //   fBeams[beam].SetOscillationParameters(0., 0., 0.);
  //   fBeams[beam].MakeResultPlots();
  //   return;

  // outer loop - global minima:

  float chisq = 0.;
  float chisqmin = 1.e12;
  float dmmin = 0.;
  float stmin = 0.;
  float fsmin = 0.;
  float normmin = 0.;
  float normfact = 0.;

  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;
      }
    }
  }

  // dmsq loop
  float dmsq = 0.;
  float s2t = 0.;
  float fster = 0.;
  float chisqminnorm = 9999.;
  float normtemp = 0.;
  float obs = 0.;
  float exp = 0.;
  int nmin = 0;
  int nmax = 0;

  TStopwatch sw;
  sw.Start();

  NCEnergyBin *bin = 0;

  for(int dm = 0; dm < NCType::kNumDeltaMSqrBins; dm++){

    dmsq = (1.*dm + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;
    dmsq *= 1.e-3;

    MSG("NCExtrapolationDP", Msg::kDebug) << "dmsq = " << dmsq << endl;

    // sin2theta loop
    for(int st = 0; st < NCType::kNumUMu3SqrBins; st++){

      s2t = (1.*st + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;
      MAXMSG("NCExtrapolationDP", Msg::kDebug, kNumUMu3SqrBins) << "s2t = "
                                                               << s2t << endl;

      // fsterile loop
      for (int fs = 0; fs < NCType::kNumUS3SqrBins; fs++){

        fster = (1.*fs + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;
        MAXMSG("NCExtrapolationDP", Msg::kDebug, kNumUS3SqrBins) << "fster = "
                                                                   << fster << endl;

//      fBeams[beam].SetOscillationParameters(s2t, dmsq, fster);

        // do chisq fit
        chisqminnorm = 9999.;
        normtemp = 0.;

        nmin = 50;
        nmax = 51;

        if(fsys == 1){
          nmin = 0;
          nmax = 101;
        }

        // *** normalisation loop ***
        for(int n = nmin; n < nmax; n++){

          normfact = 0.95+0.001*n;

          chisq = 0;

          //first loop is over NC/CC type
          for(int j = 0; j < 2; j++){
            for(int i = 0; i < NCType::kNumEnergyBinsFar; i++){

              bin = fBeams[beam].GetEnergyBin(i,j);
              obs = bin->GetSignal();
              exp = bin->GetTotalMC()*normfact;

              if(exp > 0){

                //calculate chisq
                chisq += 2*(exp-obs);
                if(obs > 0) chisq += 2*obs*log(obs/exp);

              }//end if expected value > 0

//            if(dm == 118 && fs == 27 && st == 98){
//              MSG("NCExtrapolationDP", Msg::kInfo)
//                << dmsq << " " << fster << " " << s2t << " "
//                << bin->GetBinCentralValue()
//                << " " << obs << " "
//                << exp << "/" << bin->GetMC() << " " << bin->GetMCBackground()
//                << " " << bin->GetMCTau() << " " << chisq << endl;
//            }


            }//end loop over energy bins
          }//end loop over nc/cc

          // add normalisation penalty term to chisq with 4% error
          chisq += TMath::Power((normfact-1)/0.04, 2);

          // remember minimum chisq in normalisation loop
          if(chisq < chisqminnorm){
            chisqminnorm = chisq;
            normtemp = normfact;
          }

        }

        //get the chisqr for this point in parameter space
        chiSqrPnts[dm][fs][st] = chisqminnorm;

        // remember global chisq minimum

        if(chisqminnorm < chisqmin){
          chisqmin = chisqminnorm;
          dmmin = dmsq;
          stmin = s2t;
          fsmin = fster;
          normmin = normtemp;
        }//end if lower chi^2
      } // end of fsterile loop
    } // end of sin2theta loop
  } // end of dmsq loop

  sw.Stop();
  MSG("NCExtrapolationDP", Msg::kInfo) << "grid search took "
                                       << sw.CpuTime() << " cpu, "
                                       << sw.RealTime() << " real"
                                       << endl;

  //set the oscillation parameters to the best fit values
  BestUMu3Sqr() = stmin;
  BestUS3Sqr() = fsmin;
  BestDeltaMSqr() = dmmin;
  fMinChiSqr = chisqmin;

  //gotta get the results ready.
//   fBeams[beam].SetOscillationParameters(fBestUMu3Sqr, fBestDeltaMSqr, fBestUS3Sqr);

  MSG("NCExtrapolationDP", Msg::kInfo) << "Best Sterile Fit at :"
                                       << " delta m^2 = " << dmmin
                                       << " f_s = " << fsmin
                                       << " sin^2(2theta) = " << stmin
                                       << " normalization = " << normmin
                                       << " chi^2 = " << chisqmin
                                       << endl;

  NCExtrapolation::FillDeltaChiSqrMaps(chiSqrPnts);

  return;
}

//----------------------------------------------------------------------
int NCExtrapolationDP::FindRecoEnergyBin(float enu, int beam)
{
  if(enu > NCType::kMaxEnergy) return NCType::kNumEnergyBinsNear-1;

  int bin = fBeams[beam].GetMCHistogram(NCType::kCC)->FindBin(enu) - 1;

  return bin;
}

//----------------------------------------------------------------------
//this method does the fit in the far detector
void NCExtrapolationDP::WriteResources()
{

  //fill the histograms in the NCBeam objects in fBeams vector
  //before calling WriteResources in the base class

  NCExtrapolation::WriteResources();

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* fileDir = gDirectory;
  TDirectory* saveDir = gDirectory;
  fileDir->cd("plots"+NCType::kExtrapolationNames[fExtrapolationType]);

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

  fnom->Write();
  fmaqewgt->Write();
  fmareswgt->Write();
  fdisfactwgt->Write();

  fnomcc->Write();
  fmaqewgtcc->Write();
  fmareswgtcc->Write();
  fdisfactwgtcc->Write();

  fnomnc->Write();
  fmaqewgtnc->Write();
  fmareswgtnc->Write();
  fdisfactwgtnc->Write();

  fdata->Write();
  fdatanc->Write();

  fmcnomcc->Write();
  fmcbfcc->Write();
  fmtempcc->Write();

  fmcnomnc->Write();
  fmcbfnc->Write();
  fmtempnc->Write();

  saveDir->cd();

  return;
}

---