NCBeam.cxx

Go to the documentation of this file.

//$Id: NCBeam.cxx,v 1.73 2009/12/18 01:30:21 rodriges Exp $
//
//NCBeam.cxx
//
//B. Rebel 3/2005

#include "NCUtils/Extrapolation/NCBeam.h"

#include "MessageService/MsgService.h"
#include "Conventions/SimFlag.h"
#include "TMath.h"

#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpAnalysisInfo.h"

#include "NCUtils/Extrapolation/NCEnergyBin.h"
#include "NCUtils/Extrapolation/NCCoordinateConverter.h" // SystPars
#include "NCUtils/NCOscProb.h"

#include <cassert>

using NCType::kNC;
using NCType::kCC;


ClassImp(NCBeam)

CVSID("$Id: NCBeam.cxx,v 1.73 2009/12/18 01:30:21 rodriges Exp $");

//-------------------------------------------------------------------
NCBeam::NCBeam() : fInfo(BeamType::BeamType_t(0), NC::RunUtil::ERunType(0))
{
  InitNulls();
}


//-------------------------------------------------------------------
NCBeam::NCBeam(Detector::Detector_t detector,
               NCBeam::Info beamInfo,
               bool useCC, bool useNC,
               bool fillBins) :
  fInfo(beamInfo),
  fDetector(detector),
  fUseCC(useCC),
  fUseNC(useNC),
  fFillBins(fillBins)
{
  Init();
}


//-------------------------------------------------------------------
void NCBeam::InitNulls()
{
  for(int nccc = kNC; nccc <= kCC; ++nccc){
    fDefaultMCHistogram[nccc] = 0;
    fMCHistogram[nccc] = 0;
    fMCTrueHistogram[nccc] = 0;
    fMCFitHistogram[nccc] = 0;
    fMCBackgroundHistogram[nccc] = 0;
    fDataGraph[nccc] = 0;
    fMCBackgroundHistogram[nccc] = 0;
    fMCNoFitNuMuToNuTau[nccc] = 0;
    fMCFitNuMuToNuTau[nccc] = 0;
    fMCNoFitNuMuToNuE[nccc] = 0;
    fMCFitNuMuToNuE[nccc] = 0;
    fMCNoFitNuEToNuE[nccc] = 0;
    fMCFitNuEToNuE[nccc] = 0;
    fDataHistogram[nccc] = 0;
    fMCNoFitBackgroundHistogram[nccc] = 0;
    fMCResidualHistogram[nccc] = 0;
  }
}


//-------------------------------------------------------------------
void NCBeam::Init()
{
  MSG("NCBeam", Msg::kDebug) << "NCBeam::Constructor" << endl;

  fResultHistogramsFilled = false;

  InitNulls();

  //make the energy bin objects for cc and nc events
  double binCentralValue = 0.;
  double binWidth = kEnergyBinWidthNear; //0.5;
  fNumEnergyBins = kNumEnergyBinsNear; //200;
  if(fDetector == Detector::kFar){
    fNumEnergyBins = kNumEnergyBinsFar; //23;
  }

  for(int i = 0; i < fNumEnergyBins; ++i){

    fEnergyBinBoundaries[i] = kEnergyBinsNear[i];
    binCentralValue = binWidth*(i*1. + 0.5);

    if(fDetector == Detector::kFar){
      fEnergyBinBoundaries[i] = kEnergyBinsFar[i];
      binWidth = kEnergyBinsFar[i+1] - kEnergyBinsFar[i];
      binCentralValue = binWidth*(0.5) + fEnergyBinBoundaries[i];
    }

    MSG("NCBeam", Msg::kDebug) << "energy bin " << i << " " << binCentralValue
                              << " " << binWidth << " " << fEnergyBinBoundaries[i]
                              << " " << fNumEnergyBins << endl;

    if(fUseCC || fDetector == Detector::kNear)
      fEnergyBins[kCC].push_back(new NCEnergyBin(binCentralValue, binWidth,
                                                 kCC));
    if(fUseNC || fDetector == Detector::kNear)
      fEnergyBins[kNC].push_back(new NCEnergyBin(binCentralValue, binWidth,
                                                 kNC));
  } // end for i

  fEnergyBinBoundaries[fNumEnergyBins] = kMaxEnergy;

  //make the fit result histograms
  TString name = "monteCarlo";
  TString nc = "NC";
  TString cc = "CC";
  TString fit = "Fit";
  TString residual = "Residual";
  TString noFit = "Nominal";
  TString bg = "Background";
  TString y = "Y";
  TString truth = "true";

  TString type = Detector::AsString(fDetector);
  type += GetInfo().GetDescription();

  for(int nccc = kNC; nccc <= kCC; ++nccc){
    TString ncccStr = nccc == kNC ? nc : cc;

    fMCHistogram[nccc] =
      NewVisEnergySpectrum(name+type+noFit+ncccStr, kRed);

    fMCTrueHistogram[nccc] =
      NewVisEnergySpectrum(name+type+truth+ncccStr, kBlack);
    fMCTrueHistogram[nccc]->SetXTitle("E_{Truth} (GeV)");

    fDefaultMCHistogram[nccc] =
      NewVisEnergySpectrum(name+"Default"+type+noFit+ncccStr, kBlack);

    fDataHistogram[nccc] =
      NewVisEnergySpectrum("Data_"+type+ncccStr, kBlack);
    fDataHistogram[nccc]->SetMarkerStyle(kFullSquare);

    fMCFitHistogram[nccc] =
      NewVisEnergySpectrum(name+type+fit+ncccStr, kBlue);

    fMCResidualHistogram[nccc] =
      NewVisEnergySpectrum(name+type+residual+ncccStr, kBlue);
    fMCResidualHistogram[nccc]->SetYTitle("Residual/GeV");
    fMCResidualHistogram[nccc]->SetLineStyle(9);

    fMCBackgroundHistogram[nccc] =
      NewVisEnergySpectrum(name+type+bg+ncccStr, kBlack, kBlack);

    fMCNoFitBackgroundHistogram[nccc] =
      NewVisEnergySpectrum(name+type+noFit+bg+ncccStr, kBlack, kBlack);

    fMCFitNuMuToNuTau[nccc] =
      NewVisEnergySpectrum(name+type+"NuMuToNuTauFit"+ncccStr, kBlack, kMagenta+2);

    fMCNoFitNuMuToNuTau[nccc] =
      NewVisEnergySpectrum(name+type+"NuMuToNuTauNoFit"+ncccStr, kBlack, kMagenta+2);

    fMCFitNuMuToNuE[nccc] =
      NewVisEnergySpectrum(name+type+"NuMuToNuEFit"+ncccStr, kBlack, kOrange-3);

    fMCNoFitNuMuToNuE[nccc] =
      NewVisEnergySpectrum(name+type+"NuMuToNuENoFit"+ncccStr, kBlack, kOrange-3);

    fMCFitNuEToNuE[nccc] =
      NewVisEnergySpectrum(name+type+"NuEToNuEFit"+ncccStr, kBlack, kGreen+2);

    fMCNoFitNuEToNuE[nccc] =
      NewVisEnergySpectrum(name+type+"NuEToNuENoFit"+ncccStr, kBlack, kGreen+2);

    TString data = "data";
    fDataGraph[nccc] = new TGraphAsymmErrors(fNumEnergyBins);
    fDataGraph[nccc]->SetName(data+type+ncccStr);
    fDataGraph[nccc]->SetTitle("");
    fDataGraph[nccc]->SetMarkerStyle(kFullCircle);
  } // end for nccc
}

//----------------------------------------------------------------------
NCBeam::~NCBeam()
{
  MSG("NCBeam", Msg::kDebug) << "NCBeam::Destructor" << endl;

  for(int nccc = kNC; nccc <= kCC; ++nccc){
    delete fDefaultMCHistogram[nccc];
    delete fMCHistogram[nccc];
    delete fMCTrueHistogram[nccc];
    delete fMCFitHistogram[nccc];
    delete fMCResidualHistogram[nccc];
    delete fMCBackgroundHistogram[nccc];
    delete fMCNoFitNuMuToNuTau[nccc];
    delete fMCFitNuMuToNuTau[nccc];
    delete fMCNoFitNuMuToNuE[nccc];
    delete fMCFitNuMuToNuE[nccc];
    delete fMCNoFitNuEToNuE[nccc];
    delete fMCFitNuEToNuE[nccc];
    delete fDataHistogram[nccc];
    delete fMCNoFitBackgroundHistogram[nccc];
    delete fDataGraph[nccc];

    for(unsigned int n = 0; n < fEnergyBins[nccc].size(); ++n)
      delete fEnergyBins[nccc][n];
  }
}

ClassImp(NCBeam::Info)

//----------------------------------------------------------------------
NCBeam::Info::Info(){}

//----------------------------------------------------------------------
NCBeam::Info::Info(BeamType::BeamType_t bt, NC::RunUtil::ERunType rt)
  : runType(rt)
{
  beamType = bt;

  switch(beamType){
  case BeamType::kL010z185i_i124:
  case BeamType::kL010z185i_i191:
  case BeamType::kL010z185i_i213:
  case BeamType::kL010z185i_i224:
  case BeamType::kL010z185i_i232:
  case BeamType::kL010z185i_i243:
  case BeamType::kL010z185i_i257:
  case BeamType::kL010z185i_i282:
  case BeamType::kL010z185i_i303:
  case BeamType::kL010z185i_i324:
    beamType = BeamType::kL010z185i;
    break;
  case BeamType::kL010z000i_i209:
  case BeamType::kL010z000i_i225:
  case BeamType::kL010z000i_i232:
  case BeamType::kL010z000i_i259:
  case BeamType::kL010z000i_i300:
  case BeamType::kL010z000i_i317:
  case BeamType::kL010z000i_i326:
  case BeamType::kL010z000i_i380:
    beamType = BeamType::kL010z000i;
    break;
  default:
    break;
  }
}

NCBeam::Info::~Info()
{
}

//----------------------------------------------------------------------
bool NCBeam::Info::operator==(const Info& i) const
{
  if(beamType != i.beamType) return false;
  if(runType  != i.runType)  return false;

  const unsigned int N = adjustedPars.size();
  const unsigned int M = i.adjustedPars.size();

  if(N != M) return false;

  for(unsigned int n = 0; n < N; ++n){
    if(adjustedPars[n] != i.adjustedPars[n]) return false;
    // TODO - comparing doubles, dodgy...
    const double eps = 1e-6;
    if(TMath::Abs(parShifts[n]-i.parShifts[n]) > eps) return false;
  }

  return true;
}

//----------------------------------------------------------------------
bool NCBeam::Info::operator<(const Info& i) const
{
#define COMPARE_BEAMINFO(txt) \
if(txt > i.txt) return false; \
if(txt < i.txt) return true;

  COMPARE_BEAMINFO(beamType);
  COMPARE_BEAMINFO(runType);

  const int N = adjustedPars.size();
  const int M = i.adjustedPars.size();

  if(N < M) return true;
  if(N > M) return false;

  for(int n = 0; n < N; ++n){
    COMPARE_BEAMINFO(adjustedPars[n]);

    // TODO - comparing doubles, dodgy...
    const double eps = 1e-6;
    if(TMath::Abs(parShifts[n]-i.parShifts[n]) < eps) continue;
    COMPARE_BEAMINFO(parShifts[n]);
  }

  // Identical, therefore not less-than
  return false;

#undef COMPARE_BEAMINFO
}

//----------------------------------------------------------------------
void NCBeam::Info::SetSystShifts(bool* adj, double* shift)
{
  adjustedPars.clear();
  parShifts.clear();
  for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
    if(adj[n]){
      adjustedPars.push_back(NCType::EFitParam(n));
      parShifts.push_back(shift[n]);
    }
  }
}

//----------------------------------------------------------------------
void NCBeam::Info::SetSystShifts(const vector<NCType::EFitParam>& adj,
                                 const vector<double>& shift)
{
  adjustedPars = adj;
  parShifts = shift;
}

//----------------------------------------------------------------------
void NCBeam::Info::GetSystShifts(vector<NCType::EFitParam>&  adj,
                                 vector<double>& shift) const
{
  adj = adjustedPars;
  shift = parShifts;
}

//----------------------------------------------------------------------
TString NCBeam::Info::GetDescription() const
{
  TString ret = BeamType::AsString(beamType);
  ret += " ";
  switch(runType){
  case NC::RunUtil::kRunAll: ret += "All Runs"; break;
  case NC::RunUtil::kRunI: ret += "Run I"; break;
  case NC::RunUtil::kRunII: ret += "Run II"; break;
  case NC::RunUtil::kRunIII: ret += "Run III"; break;
  default: ret += "UNKNOWN RUN!";
  }

  for(unsigned int n = 0; n < adjustedPars.size(); ++n){
    ret += " ";
    ret += NCType::kParams[adjustedPars[n]].latexName;
    ret += " = ";
    ret += TString::Format("%lf", parShifts[n]);
  }

  return ret;
}

//----------------------------------------------------------------------
bool NCBeam::Info::IsNominal() const
{
  return adjustedPars.size() == 0;
}

//----------------------------------------------------------------------
bool NCBeam::Info::IsShifted(NCType::EFitParam param, double* val) const
{
  for(unsigned int n = 0; n < adjustedPars.size(); ++n){
    if(adjustedPars[n] == param){
      if(val) *val = parShifts[n];
      return true;
    }
  }
  return false;
}

//----------------------------------------------------------------------
NC::RunUtil::ERunType NCBeam::Info::GetRunType() const
{
  return runType;
}

//----------------------------------------------------------------------
ostream& operator<<(ostream& os, const NCBeam::Info& inf)
{
  os << inf.GetDescription();
  return os;
}


//----------------------------------------------------------------------
NCEnergyBin* NCBeam::GetEnergyBin(int i, NCType::EEventType nccc) const
{
  if(int(fEnergyBins[nccc].size()) > i)
    return fEnergyBins[nccc][i];
  else
    MSG("NCBeam", Msg::kWarning) << "cant find requested bin "
                                 << "object - returning null "
                                 << "pointer" << endl;

  return 0;
}

//----------------------------------------------------------------------
//returns maximum energy for events used in analysis
double NCBeam::GetMaximumEnergy() const
{
  return kMaxEnergy;
}

//----------------------------------------------------------------------
int NCBeam::GetNumberEnergyBins(NCType::EEventType nccc) const
{
  return int(fEnergyBins[nccc].size());
}

//----------------------------------------------------------------------
void NCBeam::GetEnergyBinBoundaries(double* bins) const
{
  for(int i = 0; i < fNumEnergyBins; ++i)
    bins[i] = fEnergyBinBoundaries[i];
}

//----------------------------------------------------------------------
Detector::Detector_t NCBeam::GetDetector() const
{
  return fDetector;
}

//----------------------------------------------------------------------
BeamType::BeamType_t NCBeam::GetBeamType() const
{
  return fInfo.beamType;
}

//----------------------------------------------------------------------
NC::RunUtil::ERunType NCBeam::GetRunType() const
{
  return fInfo.GetRunType();
}

//----------------------------------------------------------------------
NCBeam::Info NCBeam::GetInfo() const
{
  return fInfo;
}

//----------------------------------------------------------------------
void NCBeam::SetRunType(NC::RunUtil::ERunType runType)
{
  fInfo.runType = runType;
}

//----------------------------------------------------------------------
//if you want to use MC as data then just change the recoInfo->weight
//value appropriately before passing it into this method - works for
//changes in systematic parameters as well as oscillations
void NCBeam::AddEvent(ANtpHeaderInfo* headerInfo,
                      ANtpBeamInfo* /*beamInfo*/,
                      ANtpRecoInfo* recoInfo,
                      ANtpAnalysisInfo* analysisInfo,
                      ANtpTruthInfoBeam* truthInfo,
                      bool useMCAsData,
                      NCType::EFileType fileType)
{

  MSG("NCBeam", Msg::kDebug) << "In NCBeam::AddEvent" << endl;

  int mcCC = 0, mcNC = 0;

  //want to fill both nc and cc events in the near detector, but not necessarily
  //in the far detector
  const bool useNC = (fUseNC || fDetector == Detector::kNear);
  const bool useCC = (fUseCC || fDetector == Detector::kNear);

  //really no need to check on snarl or event quality as those events that don't pass
  //the cuts dont get put into the ntuples from the extraction code.

  //which reconstructed energy bin is the event in?
  int bin = FindEnergyBin(recoInfo->nuEnergy);

  if(bin < fNumEnergyBins && bin > -1){
    if(headerInfo->dataType == SimFlag::kData || (truthInfo && useMCAsData)){
//       MAXMSG("NCBeam", Msg::kInfo, 20) << "DATA " << headerInfo->detector << "/"
//                                     << (int)Detector::kNear << " "
//                                     << headerInfo->dataType << "/"
//                                     << (int)SimFlag::kData << " "
//                                     << analysisInfo->isNC << " "
//                                     << analysisInfo->isCC << " "
//                                     << headerInfo->run << " "
//                                     << headerInfo->snarl << " "
//                                     << endl;
      if(analysisInfo->isNC > 0 && useNC){
        ++mcNC;
        if(fFillBins) fEnergyBins[kNC][bin]->AddEventToBin(recoInfo);
      }
      else if(analysisInfo->isCC > 0 && useCC){
        ++mcCC;
        if(fFillBins) fEnergyBins[kCC][bin]->AddEventToBin(recoInfo);
      }
    }//end if data
    else if(headerInfo->dataType == SimFlag::kMC){

//       if(fileType == NCType::kElectronFile){
//      MAXMSG("NCBeam", Msg::kInfo, 20)
//        << "MC " << headerInfo->run << " " << headerInfo->snarl << " "
//        << truthInfo->interactionType << " "
//        << recoInfo->weight << " "
//        << energy << " " << recoInfo->nuEnergy << " "
//        << FindEnergyBin(energy) << endl;
//       }

      if(analysisInfo->isNC > 0 && useNC){
        ++mcNC;
        if(fFillBins) fEnergyBins[kNC][bin]->AddEventToBin(recoInfo, truthInfo, fileType);
        if(truthInfo->nuFlavor == 14)
          fDefaultMCHistogram[kNC]->Fill(recoInfo->nuEnergy, recoInfo->weight);
      }
      else if(analysisInfo->isCC > 0 && useCC){
        ++mcCC;
        if(fFillBins) fEnergyBins[kCC][bin]->AddEventToBin(recoInfo, truthInfo, fileType);
        if(truthInfo->nuFlavor == 14)
          fDefaultMCHistogram[kCC]->Fill(recoInfo->nuEnergy, recoInfo->weight);
      }

    }//end if monte carlo
  }//end if acceptable value for energy bin

  MSG("NCBeam", Msg::kDebug) << fInfo.beamType << " mc cc = " << mcCC
                             << " nc = " << mcNC << endl;
}

//----------------------------------------------------------------------
int NCBeam::FindEnergyBin(double energy)
{
  //do a binary search to find the appropriate bin - the first entry in binBoundaries
  //is 0. as the binary search returns the nearest element smaller than the
  //desired value
  if(energy > kMaxEnergy) return fNumEnergyBins;

  return TMath::BinarySearch(fNumEnergyBins, fEnergyBinBoundaries, energy);
}

//----------------------------------------------------------------------
void NCBeam::MakeResultPlots()
{
//   MSG("NCBeam", Msg::kInfo) << "NCBeam::MakeResultPlots() "
//                          << " " << Detector::AsString(fDetector)
//                          << " " << BeamType::AsString(fInfo.beamType)
//                          << " " << fInfo.runType << endl;

  for(int nccc = kNC; nccc <= kCC; ++nccc){
    MakeMCHistogram(fMCHistogram[nccc],
                    fMCTrueHistogram[nccc],
                    fMCNoFitBackgroundHistogram[nccc],
                    fMCNoFitNuMuToNuTau[nccc],
                    fMCNoFitNuMuToNuE[nccc],
                    fMCNoFitNuEToNuE[nccc],
                    fMCFitHistogram[nccc],
                    fMCBackgroundHistogram[nccc],
                    fMCFitNuMuToNuTau[nccc],
                    fMCFitNuMuToNuE[nccc],
                    fMCFitNuEToNuE[nccc]);

    MakeResidualHistogram(fEnergyBins[nccc], fMCFitHistogram[nccc],
                          fDataHistogram[nccc], fMCResidualHistogram[nccc]);
  }

  MakeDataGraph();
}

//----------------------------------------------------------------------
void NCBeam::MakeMCHistogram(TH1F* mc,
                             TH1F* mcTrue,
                             TH1F* noFitBg,
                             TH1F* noFitNuMuToNuTau,
                             TH1F* noFitNuMuToNuE,
                             TH1F* noFitNuEToNuE,
                             TH1F* fit,
                             TH1F* fitBg,
                             TH1F* fitNuMuToNuTau,
                             TH1F* fitNuMuToNuE,
                             TH1F* fitNuEToNuE)
{
  //MSG("NCBeam", Msg::kInfo) << "   NCBeam::MakeMCHistogram()" << endl;

  //make a histogram to hold the baseline monte carlo

  double normalization = 1.;
  TString yTitle = "Events/ GeV";
  double width = 1.;

  //for the near detector scale the events down by 1.e3 so the axis
  //label looks reasonable
  if(fDetector == (int)Detector::kNear){
    normalization = 1.e3;
    yTitle = "10^{3} Events/GeV";
  }

  if(mc->Integral() > 0. || fit->Integral() > 0.){
//     MSG("NCBeam", Msg::kInfo) << "nominal or fit histograms already filled"
//                            << " do scaling" << endl;
    //assume if the mc histogram has some entries it is
    //already filled and just needs to be scaled
    if(mc->Integral() > 0.){
      mc->Scale(1./normalization);
      mcTrue->Scale(1./normalization);
      noFitBg->Scale(1./normalization);

      mc->SetYTitle(yTitle);
      mcTrue->SetYTitle(yTitle);
      noFitBg->SetYTitle(yTitle);

      //loop over the bins and divide each by the width
      //to put the scale in terms of events/GeV
      double events = 0.;
      for(int i = 0; i < mc->GetNbinsX(); ++i){
        events = mc->GetBinContent(i+1);
        width = mc->GetBinWidth(i+1);
        mc->SetBinContent(i+1, events/width);

        MSG("NCBeam", Msg::kDebug) << "mc histogram bin " << i+1 << "/" << mc->GetNbinsX()
                                   << " events " << events << " width " << width
                                   << " " << mc->GetBinContent(i+1) << endl;
      }

      for(int i = 0; i < mcTrue->GetNbinsX(); ++i){
        events = mcTrue->GetBinContent(i+1);
        width = mcTrue->GetBinWidth(i+1);
        mcTrue->SetBinContent(i+1, events/width);

        MSG("NCBeam", Msg::kDebug) << "mc true histogram bin " << i+1 << "/" << mcTrue->GetNbinsX()
                                   << " events " << events << " width " << width
                                   << " " << mcTrue->GetBinContent(i+1) << endl;
      }

      for(int i = 0; i < noFitBg->GetNbinsX(); ++i){
        events = noFitBg->GetBinContent(i+1);
        width = noFitBg->GetBinWidth(i+1);
        noFitBg->SetBinContent(i+1, events/width);
      }

      if(noFitNuMuToNuTau){
        noFitNuMuToNuTau->Scale(1./normalization);
        noFitNuMuToNuTau->SetYTitle(yTitle);

        for(int i = 0; i < noFitNuMuToNuTau->GetNbinsX(); ++i){
          events = noFitNuMuToNuTau->GetBinContent(i+1);
          width = noFitNuMuToNuTau->GetBinWidth(i+1);
          noFitNuMuToNuTau->SetBinContent(i+1, events/width);
        }
      }//end if no fit tau

      if(noFitNuMuToNuE){
        noFitNuMuToNuE->Scale(1./normalization);
        noFitNuMuToNuE->SetYTitle(yTitle);

        for(int i = 0; i < noFitNuMuToNuE->GetNbinsX(); ++i){
          events = noFitNuMuToNuE->GetBinContent(i+1);
          width = noFitNuMuToNuE->GetBinWidth(i+1);
          noFitNuMuToNuE->SetBinContent(i+1, events/width);
        }
      }//end if no fit electron

      if(noFitNuEToNuE){
        noFitNuEToNuE->Scale(1./normalization);
        noFitNuEToNuE->SetYTitle(yTitle);

        for(int i = 0; i < noFitNuEToNuE->GetNbinsX(); ++i){
          events = noFitNuEToNuE->GetBinContent(i+1);
          width = noFitNuEToNuE->GetBinWidth(i+1);
          noFitNuEToNuE->SetBinContent(i+1, events/width);
        }
      }//end if no fit electron
    }//end if mc->Integral > 0.

    //assume if the mc histogram has some entries it is
    //already filled and just needs to be scaled
    if(fit->Integral() > 0.){
      fit->Scale(1./normalization);
      fitBg->Scale(1./normalization);

      fit->SetYTitle(yTitle);
      fitBg->SetYTitle(yTitle);

      //loop over the bins and divide each by the width
      //to put the scale in terms of events/GeV
      double events = 0.;
      double width = 1.;
      for(int i = 0; i < fit->GetNbinsX(); ++i){
        events = fit->GetBinContent(i+1);
        width = fit->GetBinWidth(i+1);
        fit->SetBinContent(i+1, events/width);
      }
      for(int i = 0; i < fitBg->GetNbinsX(); ++i){
        events = fitBg->GetBinContent(i+1);
        width = fitBg->GetBinWidth(i+1);
        fitBg->SetBinContent(i+1, events/width);
      }

      if(fitNuMuToNuTau){
        fitNuMuToNuTau->Scale(1./normalization);
        fitNuMuToNuTau->SetYTitle(yTitle);

        for(int i = 0; i < fitNuMuToNuTau->GetNbinsX(); ++i){
          events = fitNuMuToNuTau->GetBinContent(i+1);
          width = fitNuMuToNuTau->GetBinWidth(i+1);
          fitNuMuToNuTau->SetBinContent(i+1, events/width);
        }
      }//end if fitNuMuToNuTau

      if(fitNuMuToNuE){
        fitNuMuToNuE->Scale(1./normalization);
        fitNuMuToNuE->SetYTitle(yTitle);

        for(int i = 0; i < fitNuMuToNuE->GetNbinsX(); ++i){
          events = fitNuMuToNuE->GetBinContent(i+1);
          width = fitNuMuToNuE->GetBinWidth(i+1);
          fitNuMuToNuE->SetBinContent(i+1, events/width);
        }
      }//end if fitNuMuToNuE

      if(fitNuEToNuE){
        fitNuEToNuE->Scale(1./normalization);
        fitNuEToNuE->SetYTitle(yTitle);

        for(int i = 0; i < fitNuEToNuE->GetNbinsX(); ++i){
          events = fitNuEToNuE->GetBinContent(i+1);
          width = fitNuEToNuE->GetBinWidth(i+1);
          fitNuEToNuE->SetBinContent(i+1, events/width);
        }
      }//end if no fit electron

    }//end if fit->Integral > 0.

  }//end if fit or nominal histograms have values
}

//----------------------------------------------------------------------
void NCBeam::MakeDataHistogram(std::vector<NCEnergyBin*>& energyBins,
                               TH1F* data)
{
  //MSG("NCBeam", Msg::kInfo) << "   NCBeam::MakeDataHistogram()" << endl;

  //make a histogram to hold the data
  double totalEvents = 0.;
  TString yTitle = "Events/GeV";

  //for the near detector scale the events down by 1.e3 so the axis
  //label looks reasonable
  if(fDetector == (int)Detector::kNear){
    yTitle = "10^{3} Events/GeV";
  }

  //assume if the data histogram has some entries it is
  //already filled and just needs to be scaled

  //scale the size of the entries by the inverse of the
  //bin width in the graph making stage
  if(data->Integral() > 0.){
    data->SetYTitle(yTitle);
    return;
  }

  double binCenter = 0.;
  for(int i = 0; i < fNumEnergyBins; ++i){
    if(energyBins.size() > 0){
      binCenter = energyBins[i]->GetBinCentralValue();
      totalEvents = energyBins[i]->GetSignal();

      //BJR - 7/18/07 - dont divide by the bin width here.
      //that happens in MakeGraphFromHistogram()

      data->Fill(binCenter, (totalEvents));
      data->SetBinError(i+1,TMath::Sqrt(totalEvents));
    }
  }

  data->SetYTitle(yTitle);
}

//----------------------------------------------------------------------
void NCBeam::MakeResidualHistogram(std::vector<NCEnergyBin*>& energyBins,
                                   TH1F* fit,
                                   TH1F* data,
                                   TH1F* resid)
{
  TString yTitle = "Residual/0.5 GeV";

  double mcFit = 0.;
  double d = 0.;
  double err = 1.;

  //for the near detector scale the events down by 1.e3 so the axis
  //label looks reasonable
  if(fDetector == (int)Detector::kNear){
    yTitle = "10^{3} Residual/0.5 GeV";
  }

  //assume if the mc histogram has some entries it is
  //already filled and just needs to be scaled
  if(resid->Integral() > 0.) return;

  for(int i = 0; i < fNumEnergyBins; ++i){
    err = 1.;
    if(energyBins.size() > 0){
      mcFit = fit->GetBinContent(i+1);
      d = data->GetBinContent(i+1);
      if(mcFit>0.) err = TMath::Sqrt(mcFit);
      resid->Fill(data->GetBinCenter(i+1),(d-mcFit)/err);
    }
  }

  resid->SetYTitle(yTitle);
}

//----------------------------------------------------------------------
void NCBeam::MakeDataGraph()
{
  //MSG("NCBeam", Msg::kInfo) << "   NCBeam::MakeDataGraph()" << endl;

  for(int nccc = kNC; nccc <= kCC; ++nccc){
    if(fDataHistogram[nccc]->Integral() < 1.)
      MakeDataHistogram(fEnergyBins[nccc], fDataHistogram[nccc]);

    MakeGraphFromHistogram(fDataHistogram[nccc], fDataGraph[nccc]);
  }
}

//----------------------------------------------------------------------
void NCBeam::MakeGraphFromHistogram(TH1F* hist,
                                    TGraphAsymmErrors* graph) const
{
  //MSG("NCBeam", Msg::kInfo) << "   NCBeam::MakeGraphFromHistogram()" << endl;

  //make a histogram to hold the baseline monte carlo
  double normalization = 1.;
  TString yTitle = "Events/GeV";

  //for the near detector scale the events down by 1.e3 so the axis
  //label looks reasonable
  if(fDetector == Detector::kNear){
    normalization = 1.e3;
    yTitle = "10^{3} Events/GeV";
  }

  const double poissonErrorsUp[] = { 1.29,  2.75,  4.25,  5.30,  6.78,
                                     7.81,  9.28, 10.30, 11.32, 12.79,
                                     13.81, 14.82, 16.29, 17.30, 18.32,
                                     19.32, 20.80, 21.81, 22.82, 23.82};
  const double poissonErrorsDown[] = { 0.00,  0.37,  0.74,  1.10,  2.34,
                                       2.75,  3.82,  4.25,  5.30,  6.33,
                                       6.78,  7.81,  8.83,  9.28, 10.30,
                                       11.32, 12.33, 12.79, 13.81, 14.82};

  for(int i = 0; i < fNumEnergyBins; ++i){
    const double norm = normalization*hist->GetBinWidth(i+1);
    const double x = hist->GetBinCenter(i+1);
    double y = hist->GetBinContent(i+1);

    MSG("NCBeam", Msg::kDebug) << hist->GetName() << " "
                               << x << " " << y << endl;

    double errorYLow, errorYHigh;

    if(y < 1.){
      errorYLow = y - poissonErrorsDown[0];
      errorYHigh = poissonErrorsUp[0] - y;
    }
    else if(y < 19.){
      errorYLow = y - poissonErrorsDown[TMath::Nint(y)];
      errorYHigh = poissonErrorsUp[TMath::Nint(y)] - y;
    }
    else{
      errorYLow = TMath::Sqrt(y);
      errorYHigh = errorYLow;
    }

    y /= norm;
    errorYLow /= norm;
    errorYHigh /= norm;

    graph->SetPoint(i, x, y);
    graph->SetPointEYhigh(i, errorYHigh);
    graph->SetPointEYlow(i, errorYLow);

    MSG("NCBeam", Msg::kDebug) << y << " " << TMath::Nint(y)
                               << " " << poissonErrorsDown[TMath::Nint(y)]
                               << " " << poissonErrorsUp[TMath::Nint(y)]
                               << " " << errorYLow << " " << errorYHigh
                               << endl;

  }//end loop over energy bins

  graph->GetHistogram()->SetXTitle("E_{vis} (GeV)");
  graph->GetHistogram()->SetYTitle(yTitle);
  graph->GetHistogram()->GetXaxis()->CenterTitle();
  graph->GetHistogram()->GetYaxis()->CenterTitle();
}

//----------------------------------------------------------------------
const TH1F* NCBeam::GetDefaultMCHistogram(NCType::EEventType nccc) const
{
  return fDefaultMCHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCHistogram(NCType::EEventType nccc)
{
  return fMCHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCTrueHistogram(NCType::EEventType nccc)
{
  return fMCTrueHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetDataHistogram(NCType::EEventType nccc)
{
  return fDataHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCFitHistogram(NCType::EEventType nccc)
{
  return fMCFitHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCResidualHistogram(NCType::EEventType nccc)
{
  return fMCResidualHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCBackgroundHistogram(NCType::EEventType nccc)
{
  return fMCBackgroundHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCNoFitBackgroundHistogram(NCType::EEventType nccc)
{
  return fMCNoFitBackgroundHistogram[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCFitNuMuToNuTauHistogram(NCType::EEventType nccc)
{
  return fMCFitNuMuToNuTau[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCNoFitNuMuToNuTauHistogram(NCType::EEventType nccc)
{
  return fMCNoFitNuMuToNuTau[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCFitNuMuToNuEHistogram(NCType::EEventType nccc)
{
  return fMCFitNuMuToNuE[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCNoFitNuMuToNuEHistogram(NCType::EEventType nccc)
{
  return fMCNoFitNuMuToNuE[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCFitNuEToNuEHistogram(NCType::EEventType nccc)
{
  return fMCFitNuEToNuE[nccc];
}

//----------------------------------------------------------------------
TH1F* NCBeam::GetMCNoFitNuEToNuEHistogram(NCType::EEventType nccc)
{
  return fMCNoFitNuEToNuE[nccc];
}

//----------------------------------------------------------------------
TGraphAsymmErrors* NCBeam::GetDataGraph(NCType::EEventType nccc)
{
  return fDataGraph[nccc];
}

//----------------------------------------------------------------------
void NCBeam::Reset(bool data, bool mc)
{
  for(int nccc = kNC; nccc <= kCC; ++nccc){
    if(data){
      fDataHistogram[nccc]->Reset("ICE");
      delete fDataGraph[nccc];
      fDataGraph[nccc] = new TGraphAsymmErrors;
    }
    if(mc){
      fDefaultMCHistogram[nccc]->Reset("ICE");
      fMCHistogram[nccc]->Reset("ICE");
      fMCTrueHistogram[nccc]->Reset("ICE");
      fMCFitHistogram[nccc]->Reset("ICE");
      fMCResidualHistogram[nccc]->Reset("ICE");
      fMCBackgroundHistogram[nccc]->Reset("ICE");
      fMCNoFitNuMuToNuTau[nccc]->Reset("ICE");
      fMCFitNuMuToNuTau[nccc]->Reset("ICE");
      fMCNoFitNuMuToNuE[nccc]->Reset("ICE");
      fMCFitNuMuToNuE[nccc]->Reset("ICE");
      fMCNoFitNuEToNuE[nccc]->Reset("ICE");
      fMCFitNuEToNuE[nccc]->Reset("ICE");
      fMCNoFitBackgroundHistogram[nccc]->Reset("ICE");
    }

    for(unsigned int n = 0; n < fEnergyBins[nccc].size(); ++n)
      fEnergyBins[nccc][n]->Reset(data, mc);
  }
}

//----------------------------------------------------------------------
//write out all graphs and histograms
void NCBeam::WriteResources()
{
  WritePredictionResources();

  for(int nccc = kNC; nccc <= kCC; ++nccc){
    fDataHistogram[nccc]->Write();
    fDataGraph[nccc]->Write();
  }
}

//----------------------------------------------------------------------
//write out all graphs and histograms
void NCBeam::WritePredictionResources()
{
  for(int nccc = kNC; nccc <= kCC; ++nccc){
    fMCHistogram[nccc]->Write();
    fMCTrueHistogram[nccc]->Write();
    fDefaultMCHistogram[nccc]->Write();
    fMCFitHistogram[nccc]->Write();
    fMCResidualHistogram[nccc]->Write();
    fMCBackgroundHistogram[nccc]->Write();
    fMCNoFitNuMuToNuTau[nccc]->Write();
    fMCFitNuMuToNuTau[nccc]->Write();
    fMCNoFitNuMuToNuE[nccc]->Write();
    fMCFitNuMuToNuE[nccc]->Write();
    fMCNoFitNuEToNuE[nccc]->Write();
    fMCFitNuEToNuE[nccc]->Write();
    fMCNoFitBackgroundHistogram[nccc]->Write();
  }
}

//----------------------------------------------------------------------
/*
void NCBeam::Add(const NCBeam& beam)
{
  const bool useNC = (fUseNC || fDetector == Detector::kNear);
  const bool useCC = (fUseCC || fDetector == Detector::kNear);

  if(beam.GetBeamType() == fInfo.beamType &&
     beam.GetDetector() == fDetector){

    fDataPOT += beam.GetDataPOT();
    fMCPOT += beam.GetMCPOT();

    //loop over the energy bins in the beam and add them to the bins in this beam
    for(unsigned int i = 0; i < fEnergyBins[kCC].size() && useCC; ++i)
      fEnergyBins[kCC][i]->Add(beam.GetEnergyBin(i, kCC));

    for(unsigned int i = 0; i < fEnergyBins[kNC].size() && useNC; ++i)
      fEnergyBins[kNC][i]->Add(beam.GetEnergyBin(i, kNC));

    // Need to also add histograms. This is my random guess as to the ones
    // that we calculate internally instead of being filled externally.
    for(int nccc = kNC; nccc <= kCC; ++nccc){
      fDefaultMCHistogram[nccc]->Add(fDefaultMCHistogram[nccc]);
      fMCHistogram[nccc]->Add(fMCHistogram[nccc]);
      fMCTrueHistogram[nccc]->Add(fMCTrueHistogram[nccc]);
      fDataHistogram[nccc]->Add(fDataHistogram[nccc]);
    // Can't add a TGraphAsymmErrors...
    //    fDataGraph[nccc]->Add(fDataGraph[nccc]);
    }
  }
  else
    MSG("NCBeam", Msg::kWarning) << " trying to add incompatible beams "
                                 << BeamType::AsString(fInfo.beamType) << "/"
                                 << BeamType::AsString(beam.GetBeamType())
                                 << " " << Detector::AsString(fDetector)
                                 << "/" << Detector::AsString(beam.GetDetector())
                                 << endl;
}

*/
//----------------------------------------------------------------------
TH1F* NCBeam::NewVisEnergySpectrum(TString name, int col, int fillcol) const
{
  TH1F* ret = new TH1F(name, "", fNumEnergyBins, fEnergyBinBoundaries);
  ret->SetXTitle("E_{vis} (GeV)");
  ret->SetYTitle("Events/GeV");
  ret->GetXaxis()->CenterTitle();
  ret->GetYaxis()->CenterTitle();
  ret->SetLineColor(col);

  if(fillcol >= 0){
    ret->SetFillColor(fillcol);
    ret->SetFillStyle(3004);
  }

  return ret;
}


//----------------------------------------------------------------------
void NCBeam::FillResultHistograms(const NC::OscProb::OscPars* oscPars,
                                  const NC::SystPars* systPars)
{
  if(fResultHistogramsFilled){
    MSG("NCBeam", Msg::kWarning) << "Calling FillResultHistograms again on "
                                 << "the same beam "
                                 << GetInfo().GetDescription()
                                 << GetDetector()
                                 << ". Call ignored." << endl;
    return;
  }
  fResultHistogramsFilled = true;

  using NCType::EEventType;
  using NCType::kNC;
  using NCType::kCC;

  // These are the variables we load the event information into
  double trueEnergy, recoShowerE, recoMuonE, weight;
  int flavor;

  const Detector::Detector_t det = GetDetector();
  const double baseLine = (det == Detector::kFar) ? NCType::kBaseLineFar
                                                  : NCType::kBaseLineNear;

  //loop over NC/CC energy bins in beam
  for(EEventType nccc = kNC; nccc <= kCC; nccc = EEventType(int(nccc)+1)){

    if(nccc == kNC && det == Detector::kFar && !fUseNC) continue;
    if(nccc == kCC && det == Detector::kFar && !fUseCC) continue;

    const EEventType sigType = EEventType(nccc);
    const EEventType bkgType = (nccc == kNC) ? kCC : kNC;

    const int I = GetNumberEnergyBins(nccc);
    for(int i = 0; i < I; ++i){
      NCEnergyBin* bin = GetEnergyBin(i, nccc);

      const int mcSize             = bin->GetMCSignalVectorSize();
      const int mcSize_bkg         = bin->GetMCBackgroundVectorSize();
      const int mcSize_NuMuToNuTau = bin->GetMCNuTauVectorSize();
      const int mcSize_NuEToNuE    = bin->GetMCBeamNuEVectorSize();
      const int mcSize_NuMuToNuE   = bin->GetMCOscNuEVectorSize();

      for(int e = 0; e < mcSize; ++e){
        bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
        double survivalProb = oscPars->TransitionProbability(NCType::kNuMuToNuMu,
                                                             sigType,
                                                             baseLine,
                                                             trueEnergy);

        GetMCHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());

        GetMCTrueHistogram(sigType)->Fill(trueEnergy,weight);

        GetMCFitHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());
      }//end loop over signal

      for(int e = 0; e < mcSize_bkg; ++e){
        bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);

        double survivalProb = oscPars->TransitionProbability(NCType::kNuMuToNuMu,
                                                             bkgType,
                                                             baseLine,
                                                             trueEnergy);

        GetMCBackgroundHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());

        GetMCHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());

        GetMCFitHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());
      }//end loop over background

      //nutau
      for(int e = 0; e < mcSize_NuMuToNuTau; ++e){
        bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);

        double survivalProb = oscPars->TransitionProbability(NCType::kNuMuToNuTau,
                                                             kCC,
                                                             baseLine,
                                                             trueEnergy);

        GetMCFitNuMuToNuTauHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());

        GetMCFitHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());

        GetMCNoFitNuMuToNuTauHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());
      }//end loop over taus

      //nue
      for(int e = 0; e < mcSize_NuEToNuE; ++e){
        bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);

        double survivalProb = oscPars->TransitionProbability(NCType::kNuEToNuE,
                                                             kCC,
                                                             baseLine,
                                                             trueEnergy);
        GetMCHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());

        GetMCFitHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale()
               *survivalProb);

        GetMCNoFitNuEToNuEHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());

        GetMCFitNuEToNuEHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale()
               *survivalProb);
      }//end loop over electrons

      //nue appearance
      for(int e = 0; e < mcSize_NuMuToNuE; ++e){
        bin->GetMCOscNuEInformation(trueEnergy, recoShowerE,
                                    recoMuonE, weight, flavor, e);

        double survivalProb = oscPars->TransitionProbability(NCType::kNuMuToNuE,
                                                             kCC,
                                                             baseLine,
                                                             trueEnergy);

        GetMCFitHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());

        GetMCNoFitNuMuToNuEHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*systPars->NormScale());

        GetMCFitNuMuToNuEHistogram(sigType)->
          Fill(recoShowerE*systPars->ShowerScale(det)
               + recoMuonE*systPars->TrackScale(),
               weight*survivalProb*systPars->NormScale());
      }//end loop over electron appearance

    }//end loop over energy bins

  }//end loop over nc/cc
}

---