PIDSpectrum.cxx

Go to the documentation of this file.

#include "NCUtils/Extrapolation/PIDSpectrum.h"

#include "NCUtils/Extrapolation/NCBeam.h" // for binning scheme

#include "TCanvas.h"
#include "TLatex.h"
#include "TMath.h"
#include "TPad.h"
#include "TStyle.h"
#include "TLegend.h"

#include <cmath>
#include <iostream>
#include <algorithm>

using namespace std;

ClassImp(PIDSpectrum)

PIDSpectrum::PIDSpectrum(string _name, string _title,
                         int _nPIDBins, double _PIDmin, double _PIDmax,
                         int _trueBinFactor, NCBeam::Info _beamInfo,
                         double _cutValue)
: TNamed(_name, _title),
  fNumPIDBins(_nPIDBins),
  fPIDmin(_PIDmin), fPIDmax(_PIDmax),
  fTrueBinFactor(_trueBinFactor),
  fBeamInfo(_beamInfo)
{
  fEvtTypes.push_back(kNC);
  fEvtTypes.push_back(kNumu);
  fEvtTypes.push_back(kBeamNue);
  fEvtTypes.push_back(kOscNue);
  fEvtTypes.push_back(kNutau);

  fDrawColors.resize(fEvtTypes.size());
  fDrawColors[kNC]=kRed;
  fDrawColors[kNumu]=kBlue;
  fDrawColors[kBeamNue]=kGreen+2;
  fDrawColors[kOscNue]=kOrange-3;
  fDrawColors[kNutau]=kMagenta+2;

  // Set up the binning for true energy axes
  Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar,
                                   kEnergyBinsFar);
  int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor;

  // Gah, have to construct this manually, because ROOT won't let me
  // mix-and-match the auto and manual binning for TH3s
  if (fNumPIDBins==-1) {
    // We're emulating the FarNear method - two bins, separated at the cut value
    fPIDBins=new Double_t[3];
    fPIDBins[0]=fPIDmin; fPIDBins[2]=fPIDmax;
    fPIDBins[1]=_cutValue;
    // HACK: we probably shouldn't change this from what was passed
    // in, but it makes the rest of the code easier
    fNumPIDBins=2;
  } else {
    fPIDBins=new Double_t[fNumPIDBins+1];
    for (int i=0; i<fNumPIDBins+1; ++i)
      fPIDBins[i]=fPIDmin+float(i)*(fPIDmax-fPIDmin)/fNumPIDBins;
  }

  // This is the one place where we *should* use this variable
  fNumBinEdgesDontUseThis=fNumPIDBins+1;

  for (vector<evtType>::iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
    string tmp("hTrue");
    tmp+=GetName();
    tmp+=EvtTypeAsString(*it);
    //cout << "Creating true energy hist" << endl;
    fCmpTrue[*it]=new TH2F(tmp.c_str(), "True energy spectrum;True Energy (GeV);PID",
                          nTrueEBins, trueEBins,
                          fNumPIDBins, fPIDBins);
    tmp="t2r_";
    tmp+=GetName();
    tmp+=EvtTypeAsString(*it);
    //cout << "Creating true to reco hist" << endl;
    fTrueToReco[*it]=new TH3F(tmp.c_str(), "True-to-reco;Reco Energy;PID;True energy",
                             kNumEnergyBinsFar, kEnergyBinsFar,
                             fNumPIDBins, fPIDBins,
                             nTrueEBins, trueEBins);
  }

  delete[] trueEBins;

  string tmp("hData");
  tmp+=GetName();
  //cout << "Creating data hist" << endl;
  fData=new TH2F(tmp.c_str(), "Far Data;Reco Energy;PID",
                 kNumEnergyBinsFar, kEnergyBinsFar,
                 fNumPIDBins, fPIDBins);

  tmp="hNearMC";
  tmp+=GetName();
  fNearMC=new TH2F(tmp.c_str(), "Near MC;Reco Energy;PID",
                   kNumEnergyBinsFar, kEnergyBinsFar,
                   fNumPIDBins, fPIDBins);

  tmp="hNearData";
  tmp+=GetName();
  fNearData=new TH2F(tmp.c_str(), "Near Data;Reco Energy;PID",
                     kNumEnergyBinsFar, kEnergyBinsFar,
                     fNumPIDBins, fPIDBins);

}

//======================================================================

Double_t* PIDSpectrum::MultiplyBins(int factor, int nBinsOrig,
                                    const Double_t* origBins) const
{
  Double_t* newBins=new Double_t[(nBinsOrig+1)*factor+2];

  for (int iOrig=0; iOrig<nBinsOrig+1; ++iOrig) {
    for (int j=0; j<factor; ++j) {
      newBins[iOrig*factor+j]=origBins[iOrig]+
        float(j)/factor*(origBins[iOrig+1]-origBins[iOrig]);
    }
  }

  return newBins;
}

//======================================================================

TH2F* PIDSpectrum::GetPredicted(const NC::OscProb::OscPars* coords,
                                const char* histname,
                                bool doFN,
                                TH1* nearSpectrum) const
{
  TH2F* pred=new TH2F(histname, "Predicted total",
                      kNumEnergyBinsFar, kEnergyBinsFar,
                      fNumPIDBins, fPIDBins);
  pred->SetDirectory(0);

  for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
    AddOnePredicted(pred, *it, coords);
  }

  // Because we set the bins manually using fArray, the number of
  // entries isn't calculated. This results in the histogram not
  // drawing properly, so set the number of entries manually to
  // something reasonable
  pred->SetEntries(int(pred->Integral()));

  // Multiply the predicted spectrum by the ND data/MC ratio
  if (doFN)
    FarNearCorrect(pred, fNearData, nearSpectrum ? nearSpectrum : fNearMC);

  return pred;

}

//======================================================================

TH2F* PIDSpectrum::GetOnePredicted(evtType t,
                                   const NC::OscProb::OscPars* coords,
                                   const char* histname) const
{
  // This spectrum
  string tmp = GetName();
  tmp+="_pred";
  tmp+=histname ? histname : "";
  tmp+=EvtTypeAsString(t);
  TH2F* onePred=new TH2F(tmp.c_str(),
                         Form("Predicted %s", EvtTypeAsString(t).c_str()),
                         kNumEnergyBinsFar, kEnergyBinsFar,
                         fNumPIDBins, fPIDBins);

  AddOnePredicted(onePred, t, coords);

  return onePred;
}

//======================================================================

void PIDSpectrum::AddOnePredicted(TH2F* onePred, evtType t, const NC::OscProb::OscPars* coords) const
{
  // Set up the binning for true energy axes
  Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar,
                                   kEnergyBinsFar);
  const int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor;

  TH1F* oscProb=GetOscWeightHist(t, coords, nTrueEBins, trueEBins);
  delete[] trueEBins;

  /*

  // This is some debugging code that might turn out to be useful again one day

  if (oscProb->Integral()==0 && (t==kNumu || t==kNC)) {
    cout << "oscprob empty for type " << EvtTypeAsString(t) << endl;
    cout << "Coords: ";
    for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
    cout << endl;
  }

  static bool doneAlready=false;

  if (t==kNC && !doneAlready) { oscProb->Write(); doneAlready=true; }

  */
  // Convolve it all together

  // Need to use find() because map::operator[] isn't const
  TH3F* t2r=fTrueToReco.find(t)->second;
  TH2F* trueCmp=fCmpTrue.find(t)->second;
  /*
  if ((t==kNC || t==kNumu) && hit->second->Integral()==0) {
    cout << "t2r empty for type " << EvtTypeAsString(t) << endl;
    cout << "Coords: ";
    for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
    cout << endl;
  }

  if ((t==kNC || t==kNumu) && hit2->second->Integral()==0) {
    cout << "fCmpTrue empty for type " << EvtTypeAsString(t) << endl;
    cout << "Coords: ";
    for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
    cout << endl;
  }
  */

  // TH2 and TH3 convert to 1D bin numbers when asked for a bin by 2D
  // number. It turns out that this is rather slow (probably since
  // it's all implemented by the TH1 base class). Since I know I won't
  // be accessing past the end of the range, and I know what dimension
  // the histogram has, I can calculate this manually for (hopefully)
  // a speedup

  // GetBinContent does bounds-checking, which we don't need,
  // since the loop limits ensure that we only access real
  // bins. Instead, use the internal array fArray, which for
  // some reason is public, to access the bin values without
  // incurring the time penalty from bounds-checking
  //
  // Similarly SetBinContent can be replaced by a direct write to fArray
  // The disadvantage here is that bin errors etc won't be set, but we
  // don't need them here

  // The arrangement of t2r in memory in fArray is my trueE, PID, recoE
  // most significant first. By putting the loops in this same order
  // we consider every point in the order they are layed out in memory.
  //
  // onePred has the same arrangement, except it doesn't have a trueE direction
  //
  // It is important we loop through every single bin (including underflow
  // and overflow) or the counting of the index will go wrong.

  // Begin at the start of the true-to-reco array
  float* pt2r = t2r->fArray;

  for(int iTrueE = 0; iTrueE <= nTrueEBins+1; ++iTrueE){
    const float oscWeight = oscProb->fArray[iTrueE];

    // For each new true energy go back to the start of the output array
    float* pOnePred = onePred->fArray;

    for(int iPID=0; iPID <= fNumPIDBins+1; ++iPID){
      // The trueCmp histogram has its axes the wrong way round, so have
      // to do this calculation instead of being able to use the same trick
      // as for t2r and onePred
      const int trueE_PIDBin = iTrueE+(nTrueEBins+2)*iPID;

      const float pidWeight = trueCmp->fArray[trueE_PIDBin];
      const float oscPidWeight = oscWeight * pidWeight;

      for(int iRecoE = 0; iRecoE <= kNumEnergyBinsFar+1; ++iRecoE){

        const float fillValue = (*pt2r)*oscPidWeight;

        *pOnePred += fillValue;

        // Advance to the next entry in each array
        ++pt2r;
        ++pOnePred;

      } // end for iRecoE
    } // end for iPID
  } // end for iTrueE
  /*
  if ((t==kNC || t==kNumu) && onePred->Integral()<1e-6) {
    cout << "predicted spectrum empty for type " << EvtTypeAsString(t) << endl;
    cout << "Coords: ";
    for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
    cout << endl;
  }
  */
  //  delete oscProb;
}

//==============================================================================

double PIDSpectrum::GetEventWeight(evtType t, double trueE,
                                   const NC::OscProb::OscPars* coords) const
{
  switch (t) {
  case PIDSpectrum::kNC:
    // NCExtrapolationModule only passes NC events from the nominal
    // files - NC events from the tau and electron files are
    // ignored. The rationale is that there's more stats in the
    // nominal files than the others.

    // This is the reason for the bizarre interaction type argument here.
    return coords->TransitionProbability(NCType::kNuMuToNuMu,
                                         NCType::kNC,
                                         NCType::kBaseLineFar,
                                         trueE);
    break;
  case PIDSpectrum::kNumu:
    return coords->TransitionProbability(NCType::kNuMuToNuMu,
                                         NCType::kCC,
                                         NCType::kBaseLineFar,
                                         trueE);
    break;
  case PIDSpectrum::kBeamNue:
    // Assume beam nue's don't oscillate.

    // TODO: Do this right
    // Apparently this is what's done in other places in NCUtils, so
    // maybe it can stay this way
    return 1.;
    break;
  case PIDSpectrum::kOscNue:
    return coords->TransitionProbability(NCType::kNuMuToNuE,
                                         NCType::kCC,
                                         NCType::kBaseLineFar,
                                         trueE);
    break;
  case PIDSpectrum::kNutau:
    return coords->TransitionProbability(NCType::kNuMuToNuTau,
                                         NCType::kCC,
                                         NCType::kBaseLineFar,
                                         trueE);
    break;
  default:
      // Something went wrong
    cerr << "GetOnePredicted got unknown evtType " << t << endl;
    assert(0);
  }
}

//==============================================================================

TH1F* PIDSpectrum::GetOscWeightHist(evtType t,
                                    const NC::OscProb::OscPars* coords,
                                    int nBins, Double_t* binEdges) const
{
  // Oscillation probabilities for this event type as a fn of energy
  static TH1F* oscProb=new TH1F("oscProb", "Oscillation probability", nBins, binEdges);

  //  oscProb->Reset();
  //  oscProb->SetDirectory(0);


  // The "<=" means we include the overflow bin
  for (int i=1; i<=oscProb->GetNbinsX()+1; ++i) {
    // Poor man's integral over the bin
    double Emin=oscProb->GetXaxis()->GetBinLowEdge(i);
    double Emax=oscProb->GetXaxis()->GetBinUpEdge(i);
    double trueE=0.5*(Emin+Emax);
    oscProb->fArray[i] = GetEventWeight(t, trueE, coords);
  }

  return oscProb;
}

//======================================================================

double PIDSpectrum::PIDWithinRange(double origPID) const
{
  if(origPID>fPIDmax) return fPIDmax-1e-3;
  if(origPID<fPIDmin) return fPIDmin+1e-3;
  return origPID;
}

//======================================================================

void PIDSpectrum::FillMC(Detector::Detector_t det,
                         const ANtpTruthInfoBeam* t,
                         const ANtpRecoInfo* /* r */,
                         const ANtpAnalysisInfo* ana,
                         double E, double scale)
{
  evtType evt(GetEvtType(t));

  double pidVal=PIDWithinRange(ana->separationParameter);

  // Cap energy at 120 GeV
  // The -0.01 is to make sure the event goes in the 20-120 GeV
  // bin and not the overflow bin above
  E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);

  if (det==Detector::kFar) {
    fCmpTrue[evt]->Fill(t->nuEnergy, pidVal, scale);
    fTrueToReco[evt]->Fill(E, pidVal, t->nuEnergy, scale);
  } else if (det==Detector::kNear) {
    fNearMC->Fill(E, pidVal, scale);
  } else {
    cout << "Unknown detector " << det << ". Bailing" << endl;
    exit(1);
  }
}

//======================================================================

void PIDSpectrum::FillData(Detector::Detector_t det,
                           const ANtpTruthInfoBeam* /* t */,
                           const ANtpRecoInfo*  /* r */,
                           const ANtpAnalysisInfo* ana,
                           double E,  double scale)
{

  // Cap energy at 120 GeV
  // The -0.01 is to make sure the event goes in the 20-120 GeV
  // bin and not the overflow bin above
  E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);

  double pidVal=PIDWithinRange(ana->separationParameter);
  if (det==Detector::kFar) fData->Fill(E, pidVal, scale);
  else if (det==Detector::kNear) fNearData->Fill(E, pidVal, scale);
  else {
    cout << "Unknown detector " << det << ". Bailing" << endl;
    exit(1);
  }
}

//======================================================================

void PIDSpectrum::NormalizeTrueToRecos()
{
  for (vector<evtType>::iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
    TH3F* t2r = fTrueToReco[*it];
    const int maxE = t2r->GetNbinsZ()+1;
    TH2F* ct = fCmpTrue[*it];
    TH3F* denom=(TH3F*)t2r->Clone("denom");
    denom->Reset();
    denom->SetDirectory(0);
    for (int iRecoE=1; iRecoE<kNumEnergyBinsFar+1; ++iRecoE) {
      for (int iPID=1; iPID<fNumPIDBins+1; ++iPID) {
        // The "<=" means we include the overflow bin
        for (int iTrueE=1; iTrueE<=maxE; ++iTrueE) {
          denom->SetBinContent(iRecoE, iPID, iTrueE,
                               ct->GetBinContent(iTrueE, iPID));
        }
      }
    }
    t2r->Divide(denom);
    delete denom;
  }
}

//======================================================================

PIDSpectrum::evtType PIDSpectrum::GetEvtType(const ANtpTruthInfoBeam* t) const
{
  if (t->interactionType==0)     return kNC;
  else if (abs(t->nuFlavor)==14) return kNumu;
  else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)==12) return kBeamNue;
  else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)!=12) return kOscNue;
  else if (abs(t->nuFlavor)==16) return kNutau;
  else {
    cerr << "Unknown particle type " << t->nuFlavor << endl;
    assert(0);
  }
}

//======================================================================

string PIDSpectrum::EvtTypeAsString(evtType t) const
{

  switch (t) {
  case kNC:
    return "NC";
  case kNumu:
    return "NuMu";
  case kBeamNue:
    return "BeamNuE";
  case kOscNue:
    return "OscNuE";
  case kNutau:
    return "NuTau";
  default:
    cerr << "EvtTypeAsString got unknown evtType " << t << endl;
    return "Unknown";
  }
}

//======================================================================

TCanvas* PIDSpectrum::DrawSpectrum(const NC::OscProb::OscPars* coords, 
                                   TH1* systratios,
                                   string suffix) const
{
  TH2F* hPred=GetPredicted(coords, "total", false);

  if (systratios) hPred->Multiply(systratios);
  string canName("can_");
  canName+=GetName();
  canName+=suffix;
  TCanvas *c = new TCanvas(canName.c_str(), GetTitle());
  c->Divide(2,2);

  c->cd(1);
  hPred->DrawCopy("colz");
  TLatex* la=new TLatex;
  la->SetNDC();
  la->DrawLatex(0.2, 0.95, "Predicted");

  c->cd(2);
  fData->DrawCopy("colz");
  la->DrawLatex(0.2, 0.95, "Data");

  c->cd(3);
  TH2F* hRatio=(TH2F*)fData->Clone("ratio");
  hRatio->Divide(hPred);
  hRatio->Draw("colz");
  la->DrawLatex(0.2, 0.95, "Ratio");

  c->cd(4);
  fCmpTrue.find(kNC)->second->DrawCopy("colz");
  la->DrawLatex(0.2, 0.95, "True NC true E");

  //  delete hPred;

  return c;
}

//======================================================================

vector<TCanvas*> PIDSpectrum::DrawEnergySlices(const NC::OscProb::OscPars* coords,
                                               TH1* systratios,
                                               string suffix) const
{
  return DrawSlices("energy", coords, systratios, suffix);
}

//======================================================================

vector<TCanvas*> PIDSpectrum::DrawPIDSlices(const NC::OscProb::OscPars* coords,
                                            TH1* systratios,
                                            string suffix) const
{
  return DrawSlices("pid", coords, systratios, suffix);
}

//======================================================================

void PIDSpectrum::ResetMC()
{
  typedef map<evtType, TH2F*>::iterator it2;
  for(it2 it = fCmpTrue.begin(); it != fCmpTrue.end(); ++it)
    it->second->Reset();

  typedef map<evtType, TH3F*>::iterator it3;
  for(it3 it = fTrueToReco.begin(); it != fTrueToReco.end(); ++it)
    it->second->Reset();
}

//======================================================================

vector<TCanvas*> PIDSpectrum::DrawSlices(string variable,
                                         const NC::OscProb::OscPars* coords,
                                         TH1* systratios,
                                         string suffix) const
{
  // Are we drawing energy slices? (If not, we're drawing PID slices)
  assert(variable == "energy" || variable == "pid");
  bool energySlice=(variable=="energy");

  vector<TCanvas*> v;
  // Canvas for the energy-integrated-over-PID plot (or vice versa)

  // ROOT is really not happy having more than one canvas with the same name
  // Keep a counter so we can keep the canvases from separate calls to this
  // function distinct.
  static int canvN=0;
  TCanvas* cTotal = new TCanvas((variable+"Total"+suffix+Form("%d", canvN)).c_str());
  // Canvas for the actual slices
  TCanvas* cSlices = new TCanvas((variable+"Slices"+suffix+Form("%d", canvN)).c_str());
  ++canvN;

  v.push_back(cTotal);   v.push_back(cSlices);

  cSlices->cd();
  int nPlots = energySlice ? fNumPIDBins : kNumEnergyBinsFar;
  int nDiv = (int)TMath::Ceil(TMath::Sqrt(nPlots));

  // Divide with more rows than columns if possible, so we get wider graphs
  if ((nDiv-1)*nDiv>=nPlots)
    cSlices->Divide(nDiv-1, nDiv);
  else if ((nDiv+1)*(nDiv-1)>=nPlots)
    cSlices->Divide(nDiv-1, nDiv+1);
  else
    cSlices->Divide(nDiv, nDiv);

  TH2F* hTotal=GetPredicted(coords, "total", false);
  if (systratios) hTotal->Multiply(systratios);
  // In principle we might be being called with parameters that don't
  // correspond to the best fit, but this'll do for now
  hTotal->SetTitle("Best fit");
  NC::OscProb::NoOscillations oscParsNoOsc;
  TH2F* hNoOsc=GetPredicted(&oscParsNoOsc, "noosc", false);
  hNoOsc->SetTitle("No Oscillations");
  // Hold the components in a map
  map<evtType, TH2F*> cmpHists;

  for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
    // Drawing everything except NC
    if (*it==kNC) continue;
    cmpHists[*it]=GetOnePredicted(*it, coords, (variable+"Slices"+suffix).c_str());
    cmpHists[*it]->SetLineColor(fDrawColors[*it]);
  }

  hTotal->SetLineColor(kRed);
  hNoOsc->SetLineColor(kRed);   hNoOsc->SetLineStyle(2);
  fData->SetLineColor(kBlack);

  // Stupid ROOT and its stupid global nonsense.
  gStyle->SetOptStat(0);

  // Draw the spectrum integrated over the other variable
  cTotal->cd();
  TH1F* noOsc1D=(TH1F*)ProjectHist(hNoOsc, variable, "NoOsc")->DrawCopy();
  noOsc1D->SetLineStyle(2);
  if (energySlice) noOsc1D->GetXaxis()->SetRangeUser(0,19.9);
  ProjectHist(hTotal, variable, "Total")->DrawCopy("same");

  for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
    // Drawing everything except NC
    if (*it==kNC) continue;
    ProjectHist(cmpHists[*it], variable,
                string("slice")+EvtTypeAsString(*it))->DrawCopy("same");
  }
  ProjectHist(fData,  variable, "Data")->DrawCopy("same e");

  ((TPad*)gPad)->BuildLegend()->SetFillColor(kWhite);
  if (!energySlice && noOsc1D->GetEntries()) ((TPad*)gPad)->SetLogy();

  // Draw the individual slices
  for (int i=1; i<nPlots+1; ++i) {
    cSlices->cd(i);

    TH1F* noOsc1DSlice=(TH1F*)ProjectHist(hNoOsc, variable,
                                          "NoOsc", i, i)->DrawCopy();
    noOsc1DSlice->SetLineStyle(2);
    if (energySlice) noOsc1DSlice->GetXaxis()->SetRangeUser(0,19.9);

    ProjectHist(hTotal, variable, "Total", i, i)->DrawCopy("same");

    for(unsigned int n = 0; n < fEvtTypes.size(); ++n){
      // Drawing everything except NC
      if (fEvtTypes[n]==kNC) continue;

      ProjectHist(cmpHists[fEvtTypes[n]], variable,
                  EvtTypeAsString(fEvtTypes[n]),
                  i, i)->DrawCopy("same");
    }
    ProjectHist(fData, variable, "Data", i, i)->DrawCopy("same e");

    if (!energySlice && noOsc1DSlice->GetEntries()) ((TPad*)gPad)->SetLogy();
  }

  return v;
}

//======================================================================

TH1D* PIDSpectrum::ProjectHist(TH2* h, string axis, string histname,
                               int lobin, int hibin) const
{
  TH1D* ret;
  histname=axis+histname;
  if(lobin) histname+=TString::Format("%d", lobin);
  else histname+="Integ";

  if (axis=="energy") ret = h->ProjectionX(histname.c_str(), lobin, hibin);
  else if (axis=="pid") ret = h->ProjectionY(histname.c_str(), lobin, hibin);
  else assert(0 && "unknown projection requested");
  ret->SetDirectory(0);
  TH1D* clone=(TH1D*)ret->Clone();
  clone->SetDirectory(0);
  return clone;
}

//======================================================================

vector<TCanvas*> PIDSpectrum::DrawNearSpectra(TH1* systratios, string suffix) const
{
  vector<TCanvas*> ret;
  TString canvName("nearSpectra "+fBeamInfo.GetDescription());
  canvName+=suffix.c_str();
  TCanvas* c = new TCanvas(canvName, "Near detector spectra");
  c->Divide(2,2);
  c->cd(1);
  TH2F* nearMCShifted=(TH2F*)fNearMC->Clone("nearMCShifted");
  if (systratios) nearMCShifted->Multiply(systratios);
  nearMCShifted->Draw("colz");
  c->cd(2);
  fNearData->Draw("colz");
  c->cd(3);
  TH2F* hRatio=(TH2F*)fNearData->Clone("hRatio");
  hRatio->Divide(nearMCShifted);
  hRatio->Draw("colz");
  ret.push_back(c);

  nearMCShifted->SetLineColor(kRed);
  fNearData->SetLineColor(kBlack);
  // 1D spectra
  TCanvas* cE=new TCanvas(canvName+"energy_1D", "Near detector 1D spectra");
  cE->cd();
  ProjectHist(nearMCShifted, "energy", "nearMC")->DrawCopy();
  ProjectHist(fNearData, "energy", "nearData")->DrawCopy("same e");
  ret.push_back(cE);

  TCanvas* cPID=new TCanvas(canvName+"pid_1D", "Near detector 1D spectra");
  cPID->cd();
  ProjectHist(nearMCShifted, "pid", "nearMC")->DrawCopy();
  ProjectHist(fNearData, "pid", "nearData")->DrawCopy("same e");
  ret.push_back(cPID);

  return ret;
}

//======================================================================

void PIDSpectrum::FarNearCorrect(TH2F* fdPred,
                                 const TH2F* ndData,
                                 const TH1* ndMC1D) const
{
  // TODO: Replace the data and MC spectra with just the ratio, which
  // we'll only need to calculate once

  // Downcast ndMC - ugh. This should be OK because that histogram came
  // from us initially.
  TH2F* ndMC = (TH2F*)ndMC1D;

  // Use the now-standard trick of accessing bin contents directly to
  // avoid unnecessary bounds-checking
  for (int i=0; i<fdPred->fN; ++i) {
    // Don't divide by zero. This should really be left to crash, but
    // I'm too lazy at the moment
    if (ndMC->fArray[i]) fdPred->fArray[i]*=ndData->fArray[i]/ndMC->fArray[i];
  }
}

---