NuMatrixSpectrum.cxx

Go to the documentation of this file.

#include <cassert>
#include <stdexcept>
#include <fstream>
#include <sstream>

#include "TCanvas.h"
#include "TPad.h"
#include "TGaxis.h"
#include "TPave.h"
#include "TF1.h"

#include "TFile.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TLatex.h"
#include "TGraphAsymmErrors.h"

#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuUtilities.h"
#include "NtupleUtils/NuFluctuator.h"

ClassImp(NuMatrixSpectrum)

CVSID("$Id: NuMatrixSpectrum.cxx,v 1.54 2010/02/05 12:32:48 bckhouse Exp $");

double NuMatrixSpectrum::compress = 2.5;

//____________________________________________________________________72
NuMatrixSpectrum::NuMatrixSpectrum() 
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  this->fSpectrum = 0;
  this->fPoT = 0.0;
}

//____________________________________________________________________72
NuMatrixSpectrum::NuMatrixSpectrum(const TH1D& spectrum)
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  this->fSpectrum = new TH1D(spectrum);
  fSpectrum->SetDirectory(0);
  this->fPoT = 0.0;
}

//____________________________________________________________________72
NuMatrixSpectrum::NuMatrixSpectrum(const TH1D& spectrum,
                                   const Double_t pot) 
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  this->fSpectrum = new TH1D(spectrum);
  fSpectrum->SetDirectory(0);
  this->fPoT = pot;
}

//____________________________________________________________________72

NuMatrixSpectrum::NuMatrixSpectrum( const string& filename,
                                    const string& histname,
                                    const string & POTHistname)
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  // Reset these - if this function fails, we don't want the user to
  // think that it passed.
  this->fSpectrum = 0;
  this->fPoT = 0.0;
  
  // Open the file
  TFile specfile(filename.c_str(), "READ");
  if (!specfile.IsOpen()) {
    ostringstream msg;
    msg << "Failed to open file " << filename;
    MSG("NuMatrixSpectrum",Msg::kError) << msg.str() << endl;
    throw runtime_error(msg.str().c_str());
  }
  
  // Grab the spectrum and POT histogram
  TH1D *spec = (TH1D*)specfile.Get(histname.c_str());
  TH1F *pot = (TH1F*)specfile.Get(POTHistname.c_str());

  // Check we actually got these, and fail if we didn't.
  if (!spec) {
    MSG("NuMatrixSpectrum",Msg::kError)
      << "Cannot read histogram " << histname << " from file "
      << filename << endl;
    return;
  }
  if (pot) {
    this->fPoT = pot->Integral();
  } else {
    MSG("NuMatrixSpectrum",Msg::kWarning)
      << "Cannot read POT histogram " << POTHistname << " from file "
      << filename << ". Defaulting to 0 POT" << endl;
    this->fPoT = 0;
  }
  
  // Now update the internal values
  this->fSpectrum = new TH1D(*spec);
  this->fSpectrum->SetDirectory(0);
}


//____________________________________________________________________72

NuMatrixSpectrum::NuMatrixSpectrum(TFile *specfile,
                                   const string& histname,
                                   const string & POTHistname)
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  // Reset these - if this function fails, we don't want the user to
  // think that it passed.
  this->fSpectrum = 0;
  this->fPoT = 0.0;
  
  // Grab the spectrum and POT histogram
  TH1D *spec = (TH1D*)specfile->Get(histname.c_str());
  TH1F *pot = (TH1F*)specfile->Get(POTHistname.c_str());
  
  // Check we actually got these, and fail if we didn't.
  if (!spec) {
    MSG("NuMatrixSpectrum",Msg::kError)
    << "Cannot read histogram " << histname << " from file "
    << specfile->GetName() << endl;
    return;
  }
  if (pot) {
    this->fPoT = pot->Integral();
  } else {
    MSG("NuMatrixSpectrum",Msg::kWarning)
    << "Cannot read histogram " << POTHistname << " from file "
    << specfile->GetName() << endl;
    this->fPoT = 0;
  }
  
  // Now update the internal values
  this->fSpectrum = new TH1D(*spec);
  this->fSpectrum->SetDirectory(0);
}
                  
//____________________________________________________________________72
NuMatrixSpectrum::NuMatrixSpectrum(const Double_t pot)
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  this->fSpectrum = 0; //new TH1D(spectrum);
  this->fPoT = pot;
}
//____________________________________________________________________72
NuMatrixSpectrum::NuMatrixSpectrum(const NuMatrixSpectrum& original)
: PoiUp(200), PoiDn(200), DoPoisson(false)
{
  // Don't copy if the original has no spectrum
  if (original.fSpectrum) {
    this->fSpectrum = new TH1D(*(original.fSpectrum));
    fSpectrum->SetDirectory(0);
  } else {
    this->fSpectrum = 0;
  }
  this->fPoT = original.fPoT;
}

//____________________________________________________________________72
NuMatrixSpectrum::~NuMatrixSpectrum()
{
  if (fSpectrum){delete fSpectrum; fSpectrum = 0;}
}

//____________________________________________________________________72
NuMatrixSpectrum& NuMatrixSpectrum
::operator = (const NuMatrixSpectrum& original)
{
  if(this == &original) return *this;

  if(fSpectrum) delete fSpectrum;

  if(original.fSpectrum){
    fSpectrum =  new TH1D(*(original.fSpectrum));
    fSpectrum->SetDirectory(0);
  }
  else{
    fSpectrum = 0;
  }
  fPoT = original.fPoT;

  return *this;
}

//____________________________________________________________________72
const char * NuMatrixSpectrum::GetName() const
{
  if (fSpectrum) return fSpectrum->GetName();
  return "Blank";
}

//____________________________________________________________________72
void NuMatrixSpectrum::SetName(const char * _name)
{
  if (fSpectrum) fSpectrum->SetName(_name);
  else MSG("NuMatrixSpectrum", Msg::kWarning) << "Failed to set name to " << _name 
    << ".  There is not yet a spectrum in this NuMatrixSpectrum." << endl;
}

//____________________________________________________________________72
void NuMatrixSpectrum::ResetSpectrum(TH1D& spectrum)
{
  if (fSpectrum) delete fSpectrum;
  fSpectrum = new TH1D(spectrum);
}

//____________________________________________________________________72
void NuMatrixSpectrum::Multiply(const TH1D* correction)
{
  fSpectrum->Multiply(correction);
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson, assuming correction " 
    << correction->GetName() << " has 0 error." << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] *= correction->GetBinContent(i+1);
      PoiDn[i] *= correction->GetBinContent(i+1);
    }
  } 
  */
}

//____________________________________________________________________72
void NuMatrixSpectrum::Multiply(TF1* correction)
{
  fSpectrum->Multiply(correction);
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson, assuming correction " 
    << correction->GetName() << " has 0 error." << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] *= correction->Eval(fSpectrum->GetBinCenter(i+1));
      PoiDn[i] *= correction->Eval(fSpectrum->GetBinCenter(i+1));
    }
  } 
  */
}

//____________________________________________________________________72
void NuMatrixSpectrum::Multiply(const TGraph* correction)
{
  for(int i=1; i<=fSpectrum->GetNbinsX(); i++) {
    fSpectrum->SetBinContent(i,fSpectrum->GetBinContent(i) *
                             correction->Eval(fSpectrum->
                                              GetBinCenter(i),0,""));
    fSpectrum->SetBinError(i,fSpectrum->GetBinError(i) *
                           correction->Eval(fSpectrum->
                                            GetBinCenter(i),0,""));
  }
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson, assuming correction " 
    << correction->GetName() << " has 0 error." << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] *= correction->Eval(fSpectrum->GetBinCenter(i+1),0,"");
      PoiDn[i] *= correction->Eval(fSpectrum->GetBinCenter(i+1),0,"");
    }
  }      
  */
}

//____________________________________________________________________72
void NuMatrixSpectrum::Multiply(const Double_t scaleFactor)
{
  fSpectrum->Scale(scaleFactor);
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson errors by " << scaleFactor << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] *= scaleFactor;
      PoiDn[i] *= scaleFactor;
    }
  }
  */
}


//____________________________________________________________________72
void NuMatrixSpectrum::Divide(const TH1D* correction, Option_t* option)
{
  fSpectrum->Divide(fSpectrum, correction, 1, 1, option);
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson, assuming correction " 
    << correction->GetName() << " has 0 error." << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] /= correction->GetBinContent(i+1);
      PoiDn[i] /= correction->GetBinContent(i+1);
    }
  }
  */
}

//____________________________________________________________________72
void NuMatrixSpectrum::Divide(const TGraph* correction)
{
  for(int i=1; i<=fSpectrum->GetNbinsX(); i++) {
    fSpectrum->SetBinContent(i,fSpectrum->GetBinContent(i) /
                             correction->Eval(fSpectrum->
                                              GetBinCenter(i),0,""));
    fSpectrum->SetBinError(i,fSpectrum->GetBinError(i) /
                           correction->Eval(fSpectrum->
                                            GetBinCenter(i),0,""));
  }
  /* Not Correct
  if (DoPoisson) {
    MSG("NuMatrixSpectrum",Msg::kInfo) << "Scaling Poisson, assuming correction " 
    << correction->GetName() << " has 0 error." << endl;
    for (int i = 0; i < fSpectrum->GetNbinsX(); i++) {
      PoiUp[i] /= correction->Eval(fSpectrum->GetBinCenter(i+1),0,"");
      PoiDn[i] /= correction->Eval(fSpectrum->GetBinCenter(i+1),0,"");
    }
  }        
  */
}

//____________________________________________________________________72
void NuMatrixSpectrum::Divide(const Double_t scaleFactor)
{
  Double_t lowCheck = 1.0e-5;
  if (scaleFactor>lowCheck){
    Multiply(1.0/scaleFactor);
  }
  else{
    MSG("NuMatrixSpectrum",Msg::kWarning)
    << "Trying to divide by a number < " << lowCheck << endl;
    assert(false);
  }
}


//____________________________________________________________________72
void NuMatrixSpectrum::MultiplyPurity(const TH1D* correction,
                                      const Double_t correctionShift)
{
  TH1D oneMinusCorrection(*correction);
  oneMinusCorrection.Reset();
  for (Int_t bin=0; bin<=oneMinusCorrection.GetNbinsX()+1; ++bin){
    oneMinusCorrection.SetBinContent(bin,1.0);
  }
  oneMinusCorrection.Add(correction,-1.0);
  TH1D background(*fSpectrum);
  background.Multiply(&oneMinusCorrection);
  background.Scale(correctionShift);

  fSpectrum->Add(&background,-1.0);
  
  if (DoPoisson) {
    DoPoisson = false;
    MSG("NuMatrixSpectrum",Msg::kWarning) << "Poisson errors cannot handle purity corrections.  Ignoring Poisson errors." << endl;
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::DividePurity(const TH1D* correction,
                                    const Double_t correctionShift)
{
  TH1D oneOverCorrMinusOne(*correction);
  oneOverCorrMinusOne.Reset();
  for (Int_t bin=0; bin<=oneOverCorrMinusOne.GetNbinsX()+1; ++bin){
    oneOverCorrMinusOne.SetBinContent(bin,1.0);
  }
  oneOverCorrMinusOne.Divide(correction);
  for (Int_t bin=0; bin<=oneOverCorrMinusOne.GetNbinsX()+1; ++bin){
    oneOverCorrMinusOne.
      SetBinContent(bin,
                    oneOverCorrMinusOne.GetBinContent(1)-1.0);
  }



  TH1D background(*fSpectrum);
  background.Multiply(&oneOverCorrMinusOne);
  background.Scale(correctionShift);

  fSpectrum->Add(&background);

  if (DoPoisson) {
    DoPoisson = false;
    MSG("NuMatrixSpectrum",Msg::kWarning) << "Poisson errors cannot handle purity corrections.  Ignoring Poisson errors." << endl;
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::RecoToTrue(const TH2D* correction)
{
  this->YToX(correction);
}

//____________________________________________________________________72
void NuMatrixSpectrum::TrueToReco(const TH2D* correction)
{
  this->XToY(correction);
}

//____________________________________________________________________72
void NuMatrixSpectrum::GenerateSystErrors() 
{
  double ebins[] = {0,     4,     8,     12,    16,    20,    30,    50, 200};
  double errbins[]    = {0.078, 0.070, 0.055, 0.060, 0.065, 0.092, 0.10,  1.0};
  double c, energy, err1, err2, errtot;
  for (int i = 1; i <= fSpectrum->GetNbinsX(); i++) {
    energy = fSpectrum->GetBinLowEdge(i);
    c = fSpectrum->GetBinContent(i);
    err1 = fSpectrum->GetBinError(i);
    
    for (int j = 0; j < 8; j++) {
      if (energy >= ebins[j] && energy < ebins[j+1]) {
        err2 = c*errbins[j];
        errtot = sqrt(err1*err1+err2*err2);
        fSpectrum->SetBinError(i, errtot);
      }
    }
  }   
}

//____________________________________________________________________72
void NuMatrixSpectrum::GenerateErrors(Double_t scale)
{
  Multiply(scale);
  double c, err1, err2, errtot;
  for (int i = 1; i <= fSpectrum->GetNbinsX(); i++) {
    c = fSpectrum->GetBinContent(i);
    err1 = fSpectrum->GetBinError(i);
    err2 = sqrt(c);
    errtot = sqrt(err1*err1+err2*err2);
    fSpectrum->SetBinError(i, errtot);
  }     
  Divide(scale);
}

//____________________________________________________________________72
void NuMatrixSpectrum::RemoveErrors()
{
  for (int i = 0; i <= fSpectrum->GetNbinsX() + 1; i++) {
    fSpectrum->SetBinError(i, 0);
  }  
}

//____________________________________________________________________72
void NuMatrixSpectrum::BinWidthNormalize()
{
  double w0, wi, w, c, e;
  w0 = fSpectrum->GetBinWidth(1);
  
  for (int i = 0; i <= fSpectrum->GetNbinsX() + 1; i++) {
    wi = fSpectrum->GetBinWidth(i);
    w = wi/w0;
    c = fSpectrum->GetBinContent(i);
    e = fSpectrum->GetBinError(i);
    fSpectrum->SetBinContent(i, c/w);
    fSpectrum->SetBinError(i, e/w);
      
    // Make sure to also scale the poisson errors we have calculated
    if (DoPoisson && i >= 1 && i <= fSpectrum->GetNbinsX()) {
      PoiUp[i-1] /= w;
      PoiDn[i-1] /= w;
    }
  }  
}


//____________________________________________________________________72
void NuMatrixSpectrum::NuBarCompressX()
{
  //cout << "Calling NuBarCompressX:" << endl
  //<< "  Make sure this step is performed after you have rebinned into" << endl
  //<< "  your desired binning scheme AND have already called BinWidthNormalize" << endl
  //<< "  since this method changes the width of the bins." << endl;
  
  // Get Current Bins
  Double_t bins_old[99];
  Double_t bins_new[99];
  GetXaxis()->GetLowEdge(bins_old);
  int Nbins = GetXaxis()->GetNbins();
  //cout << "Nbins = " << Nbins << endl;
  // Compress bins - number of bins stays the same
  //   but higher bins become half as wide
  for (int i = 0; i <= Nbins; i++) {
    if (i == Nbins) { // Overflow bin special case
      bins_new[i] = bins_new[i-1] + 150./compress; 
    }
    else if (bins_old[i] < 20.001) {
      bins_new[i] = bins_old[i];
    }
    else {
      double width = bins_old[i] - bins_old[i-1];
      bins_new[i] = bins_new[i-1] + width/compress;
    }
  }
  
  // Reassign bin contents based on bin number, 
  // not on actual bin value.
  TH1D newHist("NewHist", "", Nbins, bins_new);
  for (int i = 1; i < Nbins+1; i++) {
    newHist.SetBinContent(i, fSpectrum->GetBinContent(i));
    newHist.SetBinError(i, fSpectrum->GetBinError(i));
  }
  
  // Store our new hist as fSpectrum
  string name = fSpectrum->GetName();
  string title = fSpectrum->GetTitle();
  ResetSpectrum(newHist);
  fSpectrum->SetName(name.c_str());
  fSpectrum->SetTitle(name.c_str());
  
  // Style away the axis and labels since they are meaningless
  // fSpectrum->SetAxisColor(kWhite); // White ticks interfere with plot
  fSpectrum->GetXaxis()->SetTickLength(0);
  fSpectrum->SetLabelSize(0);
  SetRangeUser(0, bins_new[Nbins-1] - 0.005);
}


//____________________________________________________________________72
void NuMatrixSpectrum::CompressTopBins()
{
  //cout << "Calling NuBarCompressX:" << endl
  //<< "  Make sure this step is performed after you have rebinned into" << endl
  //<< "  your desired binning scheme AND have already called BinWidthNormalize" << endl
  //<< "  since this method changes the width of the bins." << endl;
  
  // Get Current Bins
  Double_t bins_old[99];
  Double_t bins_new[99];
  GetXaxis()->GetLowEdge(bins_old);
  int Nbins = GetXaxis()->GetNbins();
  //cout << "Nbins = " << Nbins << endl;
  // Compress bins - number of bins stays the same
  //   but higher bins become half as wide
  for (int i = 0; i <= Nbins-3; i++) {
    bins_new[i] = bins_old[i];
  }
  double base = bins_new[Nbins-3];
  bins_new[Nbins-2] = base+2;
  bins_new[Nbins-1] = base+4;
  bins_new[Nbins]   = base+6;
  
  
  // Reassign bin contents based on bin number, 
  // not on actual bin value.
  TH1D newHist("NewHist", "", Nbins, bins_new);
  for (int i = 1; i < Nbins+1; i++) {
    newHist.SetBinContent(i, fSpectrum->GetBinContent(i));
    newHist.SetBinError(i, fSpectrum->GetBinError(i));
  }
  
  // Store our new hist as fSpectrum
  string name = fSpectrum->GetName();
  string title = fSpectrum->GetTitle();
  ResetSpectrum(newHist);
  fSpectrum->SetName(name.c_str());
  fSpectrum->SetTitle(name.c_str());
}

//____________________________________________________________________72
TF1* NuMatrixSpectrum::GetAxisFunc()
{
    TString form = "(x<=20)*x+(x > 20)*(20+(x-20)/";
    form += compress;
    form += ")";
    TF1 * fAxis = new TF1("fAxis",form.Data(),0,50);    
    return fAxis;
}

//____________________________________________________________________72
void NuMatrixSpectrum::DrawNuBarAxes(int lsize) 
{
  gPad->Update();
  
  double aspectRatio = (double)gPad->GetCanvas()->GetWindowWidth() / (double)gPad->GetCanvas()->GetWindowHeight();
  cout << "Found aspectRatio = " << gPad->GetCanvas()->GetWindowWidth() << " / " << gPad->GetCanvas()->GetWindowHeight();
  bool wide = aspectRatio > 7./5. - 0.002;
  //bool wide = false;
  if (wide) cout << " (wide)" << endl;
  else cout << " (narrow)" << endl;

  double xmax = gPad->GetUxmax();
  double ymin = gPad->GetUymin();
  double ymax = gPad->GetUymax();
  if (gPad->GetLogy()) {
    ymin = pow(10, ymin);
    ymax = pow(10, ymax);
  }
  
  MAXMSG("NuMatrixSpectrum",Msg::kInfo,3) << "Drawing compressed axis from 0 to " 
  << xmax << " at y = " << ymin << " and y = " << ymax << endl;
  
  //cout << "compress = " << form << endl;
  TF1 * fAxis = GetAxisFunc();
  fAxis->GetName();  // Get rid of warning
    
  int div = 510;
  //if (wide) div = 510;
  //else div = 509;
  
  TGaxis *aBotFD = new TGaxis(0,ymin,xmax,ymin,"fAxis",div);
  aBotFD->SetLabelSize(0); // Drawing manually now
  aBotFD->Draw();
  
  TGaxis *aTopFD = new TGaxis(0,ymax,xmax,ymax,"fAxis",div,"-");
  aTopFD->SetLabelSize(0);
  aTopFD->Draw();

  // Draw the axis labels manually instead
  const char *xs[] = {"0", "5", "10", "15", "20", "30", "40", "50"};
  const Double_t xp[] = {0, 5, 10, 15, 20, 30, 40, 50};
  // Now draw the labels
  TLatex *p[8];
  for (int i = 0; i < 8; i++) {
    Double_t xpos;
    xpos = xp[i] <= 20 ? xp[i] : 20 + (xp[i]-20)/compress;
    //cout << "Drawing label " << xs[i] << " at " << xpos << endl;
    // Create a TLatex in a vaguely appropriate position
    p[i] = new TLatex(xpos, ymin-0.025*(ymax-ymin), xs[i]);
    p[i]->SetTextAlign(23);
    if (lsize > 0) {
      p[i]->SetTextFont(43);
      p[i]->SetTextSize(lsize);
    }
    else {
      p[i]->SetTextFont(42);
      p[i]->SetTextSize(0.06);
    }
    p[i]->Draw();
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::CutLastBin()
{
  double maxX = GetXaxis()->GetBinLowEdge(fSpectrum->GetNbinsX()) - 0.005;  
  SetRangeUser(0, maxX);
  cout << "Setting max X to " << maxX << endl;
}


//____________________________________________________________________72
void NuMatrixSpectrum::ScaleToPot(double pot)
{
  // Don't do anything unless we actually have POT
  if (fPoT < 1e-15) {
    return;
  }

  fSpectrum->Scale(pot/fPoT);
  fPoT = pot;
}

//____________________________________________________________________72
void NuMatrixSpectrum::BlessedND()
{
  static bool firstRun = true;
  if (firstRun) {
    cout << "Formatting blessed ND plot:" << endl
    << "* Binning scheme: ND Display (6) " << endl
    << "* PoT: 1e20 " << endl
    << "* Normalized to bin width" << endl
    << "* Events per 1 GeV" << endl
    << "* X axis compressed" << endl;
    firstRun = false;
  }
  RebinToScheme(NuBinningScheme::kDisplayND);
  ScaleToPot(1.e20);
  BinWidthNormalize();
  NuBarCompressX();
  fSpectrum->SetLineWidth(3);
}

//____________________________________________________________________72
void NuMatrixSpectrum::BlessedFD(bool scale)
{
  static bool firstRun = true;
  if (firstRun) {
    cout << "Formatting blessed FD plot:" << endl
         << "* Binning scheme: DisplayFD (5) " << endl;
    if (scale) cout << "* PoT: 3.2e20 " << endl;
    else       cout << "* PoT: Unchanged (" << PoT() << ")" << endl;
    cout << "* Normalized to bin width" << endl
         << "* Events per 4 GeV" << endl
         << "* X axis compressed" << endl;
    firstRun = false;
  }
  RebinToScheme(NuBinningScheme::kDisplayFD);
  if (scale) ScaleToPot(3.2e20);
  CalcPoisson();
  BinWidthNormalize();
  NuBarCompressX();
  fSpectrum->SetLineWidth(3);
}

//____________________________________________________________________72
Double_t NuMatrixSpectrum::PoiErr(double y, bool up)
{
  // Poisson Error Table
  Double_t jimeu[51];
  Double_t jimed[51];
  jimeu[0]=  1.841;  jimed[0]= 0.000;
  jimeu[1]=  3.3  ;  jimed[1]= 0.173;
  jimeu[2]=  4.638;  jimed[2]= 0.708;
  jimeu[3]=  5.918;  jimed[3]= 1.367;
  jimeu[4]=  7.163;  jimed[4]= 2.086;
  jimeu[5]=  8.382;  jimed[5]= 2.84 ;
  jimeu[6]=  9.584;  jimed[6]= 3.62 ;
  jimeu[7]=  10.77;  jimed[7]= 4.419;
  jimeu[8]=  11.95;  jimed[8]= 5.232;
  jimeu[9]=  13.11;  jimed[9]= 6.057;
  jimeu[10]= 14.27;  jimed[10]=6.891;
  jimeu[11]= 15.42;  jimed[11]=7.734;
  jimeu[12]= 16.56;  jimed[12]=8.585;
  jimeu[13]= 17.7;   jimed[13]=9.441;
  jimeu[14]= 18.83;  jimed[14]=10.3 ;
  jimeu[15]= 19.96;  jimed[15]=11.17;
  jimeu[16]= 21.08;  jimed[16]=12.04;
  jimeu[17]= 22.2 ;  jimed[17]=12.92;
  jimeu[18]= 23.32;  jimed[18]=13.8 ;
  jimeu[19]= 24.44;  jimed[19]=14.68;
  jimeu[20]= 25.55;  jimed[20]=15.57;
  jimeu[21]= 26.66;  jimed[21]=16.45;
  jimeu[22]= 27.76;  jimed[22]=17.35;
  jimeu[23]= 28.87;  jimed[23]=18.24;
  jimeu[24]= 29.97;  jimed[24]=19.14;
  jimeu[25]= 31.07;  jimed[25]=20.03;
  jimeu[26]= 32.16;  jimed[26]=20.93;
  jimeu[27]= 33.26;  jimed[27]=21.84;
  jimeu[28]= 34.35;  jimed[28]=22.74;
  jimeu[29]= 35.45;  jimed[29]=23.65;
  jimeu[30]= 36.54;  jimed[30]=24.55;
  jimeu[31]= 37.63;  jimed[31]=25.46;
  jimeu[32]= 38.72;  jimed[32]=26.37;
  jimeu[33]= 39.80;  jimed[33]=27.28;
  jimeu[34]= 40.89;  jimed[34]=28.20;
  jimeu[35]= 41.97;  jimed[35]=29.11;
  jimeu[36]= 43.06;  jimed[36]=30.03;
  jimeu[37]= 44.14;  jimed[37]=30.94;
  jimeu[38]= 45.22;  jimed[38]=31.86;
  jimeu[39]= 46.30;  jimed[39]=32.78;
  jimeu[40]= 47.32;  jimed[40]=33.70;
  jimeu[41]= 48.36;  jimed[41]=34.62;
  jimeu[42]= 49.53;  jimed[42]=35.55;
  jimeu[43]= 50.61;  jimed[43]=36.47;
  jimeu[44]= 51.68;  jimed[44]=37.39;
  jimeu[45]= 52.76;  jimed[45]=38.32;
  jimeu[46]= 53.83;  jimed[46]=39.24;
  jimeu[47]= 54.90;  jimed[47]=40.17;
  jimeu[48]= 55.98;  jimed[48]=41.10;
  jimeu[49]= 57.05;  jimed[49]=42.02;
  jimeu[50]= 58.12;  jimed[50]=42.95;
  
  Int_t n = TMath::Nint(y);
  Double_t err;
  if (up) err = jimeu[n] - n;
  else    err = n - jimed[n];
  
  return err;
}


//____________________________________________________________________72
void NuMatrixSpectrum::CalcPoisson()
{
  DoPoisson = true;
  
  const Int_t n = fSpectrum->GetNbinsX();

  Int_t y;
  for (Int_t i=1; i<=n; i++){
    y = TMath::Nint(fSpectrum->GetBinContent(i));
    PoiUp[i-1] = PoiErr(y, true);
    PoiDn[i-1] = PoiErr(y, false);
  }
}


//____________________________________________________________________72
void NuMatrixSpectrum::DrawPoisson(Option_t* opt)
{
  if (!DoPoisson) {
    MSG("NuMatrixSpectrum", Msg::kWarning) 
    << "CalcPoisson() not already called.  Poisson errors are being calculated based on current bin values."
    << endl;
    CalcPoisson();
  }
  
  TGraphAsymmErrors* g1 = new TGraphAsymmErrors(fSpectrum);
  
    
  double max = 0;
  for (Int_t i = 0; i< g1->GetN(); i++){
    Double_t x, y;
    g1->GetPoint(i, x, y);
    //cout << "i = " << i  << ", x = " << x << ", y = " << y 
    //     << ", PoiUp["<< i << "] = " << PoiUp[i] << endl;
    g1->SetPointEYlow(i, PoiDn[i]);
    g1->SetPointEYhigh(i, PoiUp[i]);
    if (y + PoiUp[i] > max) max = y+PoiUp[i];
  }
  TString option = "axis";
  option += opt;
  fSpectrum->Draw(option);
  g1->Draw("pe1same");
}


//____________________________________________________________________72
Double_t NuMatrixSpectrum::IntegralVals(double low, double high, Option_t* option)
{
  int binx1 = fSpectrum->GetXaxis()->FindFixBin(low);
  int binx2 = fSpectrum->GetXaxis()->FindFixBin(high);
  return Integral(binx1, binx2, option);
}

//____________________________________________________________________72
void NuMatrixSpectrum::ExtrapolateNDToFD(const TH2D* correction)
{
  this->XToY(correction);
}

//____________________________________________________________________72
void NuMatrixSpectrum::XToY(const TH2D* correction)
{
  const int NX = fSpectrum->GetNbinsX();

  if(correction->GetNbinsX() != NX || correction->GetNbinsY() != NX){
    MAXMSG("NuMatrixSpectrum", Msg::kError, 10)
      << "\n*** XToY: Correction has unexpected binning ***\n*** Expected "
      << NX << "x" << NX << " got "
      << correction->GetNbinsX() << "x" << correction->GetNbinsY() << " ***\n"
      << "*** Are all your spectra and matrices in consistent binnings? ***"
      << endl << endl;
  }

  //working space:
  Double_t *val = new Double_t[NX+1];
  for(int i=1;i<=NX+1;i++){
    val[i-1] = 0;
  }
  for(int i=1;i<=NX;i++){ //loop over true
    for(int j=1;j<=NX+1;j++){ //loop over reco
      val[j-1] += ( fSpectrum->GetBinContent(i) *
                    correction->GetBinContent(i,j));
    }
  }
  for(int i=1;i<=NX+1;i++) {
    fSpectrum->SetBinContent(i,val[i-1]);
  }

  delete[] val;
}

//____________________________________________________________________72
void NuMatrixSpectrum::YToX(const TH2D* correction)
{
  const int NX = fSpectrum->GetNbinsX();

  if(correction->GetNbinsY() != NX){
    MAXMSG("NuMatrixSpectrum", Msg::kError, 10)
      << "\n*** YToX: Correction has unexpected binning ***\n*** Expected "
      << NX << "x" << NX << " got "
      << correction->GetNbinsX() << "x" << correction->GetNbinsY() << " ***\n"
      << "*** Are all your spectra and matrices in consistent binnings? ***"
      << endl << endl;
  }
  

  //working space:
  Double_t *val = new Double_t[NX+1];
  for(int i=1;i<=NX+1;i++) { val[i-1] = 0; }

  const int NY = correction->GetNbinsY();
  for(int i=1;i<=NX;i++){
    for(int j=1;j<=NY;j++){
      val[j-1] += ( fSpectrum->GetBinContent(i) * 
                    correction->GetBinContent(j,i) );
    }
  }
  for(int i=1;i<=NX;i++) {
    fSpectrum->SetBinContent(i,val[i-1]);
  }

  delete[] val;
}

//____________________________________________________________________72
void NuMatrixSpectrum::Oscillate(const Double_t dm2, const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    const Double_t energy = fSpectrum->GetBinCenter(i);
    const Double_t oscProb = NuUtilities::OscillationWeight(energy, dm2, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*oscProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*oscProb);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::OscillateHybrid(const Double_t dm2, Double_t sn2)
{
  const Double_t threshold = 0.5; // The minimum spacing between oscillation minima/maxima
  Double_t minimum, previous;
  Int_t i = 2;
  minimum = (1.27 * 2 * (NuUtilities::FDDistance()/Munits::km) * dm2) / (1*TMath::Pi());
  // cout << "Initial min: " << minimum << endl;
  previous = minimum + threshold + 1;
  
  while (previous - minimum > threshold) {
    previous = minimum;
    minimum = (1.27 * 2 * (NuUtilities::FDDistance()/Munits::km) * dm2) / (i*TMath::Pi());
    ++i;
    // cout << i-1 << ": Min: " << minimum << ", Prev: " << previous << ", diff: " << previous - minimum << endl;
  }
  
  // cout << "Calculated cutoff energy: " << previous << endl;
  
  // Now min is the energy that we stop oscillating by center and start linear averaging
  for(int bin=1;bin<=fSpectrum->GetNbinsX();bin++) {
    Double_t oscProb = 0;
    
    Double_t energy = fSpectrum->GetBinCenter(bin);
    // cout << "Bin " << bin << ": ";
    // Check our cut-off point
    if (energy > previous) {
      // oscProb = 1. - (0.5 * sn2);
      oscProb = NuUtilities::OscillationWeight(energy, dm2, sn2);
     // cout << "Simple center-bin oscillation (" << energy << " = " << oscProb << "). ";
    } else {
      
      // cout << "Linear interpolation. ";
      // Use linear interpolation for the complicated parts
      const Double_t increment=0.001;
          
      if (1==bin || fSpectrum->GetNbinsX()==bin){
        // cout << "End bins. ";
        Int_t count=0;
        for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
             energy<fSpectrum->GetBinLowEdge(bin+1);
             energy+=increment){
          oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2);
          ++count;
        }
        oscProb /= count;
        //fSpectrum->SetBinContent(bin,
        //                       oscProb*fSpectrum->GetBinContent(bin));
        //fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
      }
      else{
        Double_t denominator=0.0;
        Double_t lowGrad =
          (fSpectrum->GetBinContent(bin)-fSpectrum->GetBinContent(bin-1))/
          (fSpectrum->GetBinWidth(bin)/2.0 + 
           fSpectrum->GetBinWidth(bin-1)/2.0);
        Double_t highGrad =
          (fSpectrum->GetBinContent(bin+1)-fSpectrum->GetBinContent(bin))/
          (fSpectrum->GetBinWidth(bin+1)/2.0 +
           fSpectrum->GetBinWidth(bin)/2.0);
        for (Double_t energy=fSpectrum->GetBinLowEdge(bin) + increment/2.0;
             energy<fSpectrum->GetBinLowEdge(bin+1);
             energy+=increment){
          if (energy<fSpectrum->GetBinCenter(bin)){
            oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2)*
              (fSpectrum->GetBinContent(bin-1) +
               lowGrad*(energy-fSpectrum->GetBinCenter(bin-1)));
            denominator +=
              fSpectrum->GetBinContent(bin-1) +
              lowGrad*(energy-fSpectrum->GetBinCenter(bin-1));
          }
          else {
            oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2)*
              (fSpectrum->GetBinContent(bin) +
               highGrad*(energy-fSpectrum->GetBinCenter(bin)));
            denominator +=
              fSpectrum->GetBinContent(bin) +
              highGrad*(energy-fSpectrum->GetBinCenter(bin));
          }
        }
        if (denominator>0){
          oscProb /= denominator;
        }
        else{
          oscProb = 0;
        }
      }
      
      
      
      
      
      
      
      
      
      
    }
    
    // cout << "Prob: " << oscProb << endl;
    // Now do the 'oscillation'
    fSpectrum->
      SetBinContent(bin,fSpectrum->GetBinContent(bin)*oscProb);
    fSpectrum->
      SetBinError(bin,fSpectrum->GetBinError(bin)*oscProb);
    
  }
}
//____________________________________________________________________72
void NuMatrixSpectrum::OscillateSimpleAverage(const Double_t dm2,
                                              const Double_t sn2)
{
  for (Int_t bin=1; bin<=fSpectrum->GetNbinsX(); bin++){
    Double_t oscProb=0.0;
    Double_t increment=0.001;
    Int_t count=0;
    for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
         energy<fSpectrum->GetBinLowEdge(bin+1);
         energy+=increment){
      oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2);
      ++count;
    }
    oscProb /= count;
    fSpectrum->SetBinContent(bin,
                             oscProb*fSpectrum->GetBinContent(bin));
    fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::OscillateLinearInterp(const Double_t dm2,
                                             const Double_t sn2)
{
  for (Int_t bin=1; bin<=fSpectrum->GetNbinsX(); bin++){
    Double_t oscProb=0.0;
    Double_t increment=0.001;
    if (1==bin || fSpectrum->GetNbinsX()==bin){
      Int_t count=0;
      for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
           energy<fSpectrum->GetBinLowEdge(bin+1);
           energy+=increment){
        oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2);
        ++count;
      }
      oscProb /= count;
      //fSpectrum->SetBinContent(bin,
      //                       oscProb*fSpectrum->GetBinContent(bin));
      //fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
    }
    else{
      Double_t denominator=0.0;
      Double_t lowGrad =
        (fSpectrum->GetBinContent(bin)-fSpectrum->GetBinContent(bin-1))/
        (fSpectrum->GetBinWidth(bin)/2.0 + 
         fSpectrum->GetBinWidth(bin-1)/2.0);
      Double_t highGrad =
        (fSpectrum->GetBinContent(bin+1)-fSpectrum->GetBinContent(bin))/
        (fSpectrum->GetBinWidth(bin+1)/2.0 +
         fSpectrum->GetBinWidth(bin)/2.0);
      for (Double_t energy=fSpectrum->GetBinLowEdge(bin) + increment/2.0;
           energy<fSpectrum->GetBinLowEdge(bin+1);
           energy+=increment){
        if (energy<fSpectrum->GetBinCenter(bin)){
          oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2)*
            (fSpectrum->GetBinContent(bin-1) +
             lowGrad*(energy-fSpectrum->GetBinCenter(bin-1)));
          denominator +=
            fSpectrum->GetBinContent(bin-1) +
            lowGrad*(energy-fSpectrum->GetBinCenter(bin-1));
        }
        else {
          oscProb += NuUtilities::OscillationWeight(energy, dm2, sn2)*
            (fSpectrum->GetBinContent(bin) +
             highGrad*(energy-fSpectrum->GetBinCenter(bin)));
          denominator +=
            fSpectrum->GetBinContent(bin) +
            highGrad*(energy-fSpectrum->GetBinCenter(bin));
        }
      }
      if (denominator>0){
        oscProb /= denominator;
      }
      else{
        oscProb = 0;
      }
    }
  fSpectrum->SetBinContent(bin,oscProb*fSpectrum->GetBinContent(bin));
  fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::InverseOscillate(const Double_t dm2,
                                        const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    Double_t energy = fSpectrum->GetBinCenter(i);
    Double_t oscProb = 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*oscProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*oscProb);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::InverseOscillateHybrid(const Double_t dm2, Double_t sn2)
{
  // Calculate the cut-off point
  const Double_t threshold = 0.25; // The minimum spacing between oscillation minima/maxima
  Double_t minimum, previous;
  Int_t i = 2;
  minimum = (1.27 * 2 * (NuUtilities::FDDistance()/Munits::km) * dm2) / (1*TMath::Pi());
  // cout << "Initial min: " << minimum << endl;
  previous = minimum + threshold + 1;
  
  while (previous - minimum > threshold) {
    previous = minimum;
    minimum = (1.27 * 2 * (NuUtilities::FDDistance()/Munits::km) * dm2) / (i*TMath::Pi());
    ++i;
    // cout << i-1 << ": Min: " << minimum << ", Prev: " << previous << ", diff: " << previous - minimum << endl;
  }
  
  // Now min is the energy that we stop oscillating and start normalising
  for(int bin=1;bin<=fSpectrum->GetNbinsX();bin++) {
    Double_t energy = fSpectrum->GetBinCenter(bin);
    Double_t oscProb = 0;
    
    const Double_t increment=0.001;
        
    // Check our cut-off point
    if (energy > previous) {
      oscProb = 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
    } else {
      
      
      if (1==bin || fSpectrum->GetNbinsX()==bin){
        Int_t count=0;
        for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
             energy<fSpectrum->GetBinLowEdge(bin+1);
             energy+=increment){
          oscProb += 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
          ++count;
        }
        oscProb /= count;
        //fSpectrum->SetBinContent(bin,
        //        oscProb*fSpectrum->GetBinContent(bin));
        //fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
      }
      else{
        Double_t denominator=0.0;
        Double_t lowGrad =
          (fSpectrum->GetBinContent(bin)-fSpectrum->GetBinContent(bin-1))/
          (fSpectrum->GetBinWidth(bin)/2.0 + 
           fSpectrum->GetBinWidth(bin-1)/2.0);
        Double_t highGrad =
          (fSpectrum->GetBinContent(bin+1)-fSpectrum->GetBinContent(bin))/
          (fSpectrum->GetBinWidth(bin+1)/2.0 +
           fSpectrum->GetBinWidth(bin)/2.0);
        for (Double_t energy=fSpectrum->GetBinLowEdge(bin) + increment/2.0;
             energy<fSpectrum->GetBinLowEdge(bin+1);
             energy+=increment){
          if (energy<fSpectrum->GetBinCenter(bin)){
            oscProb += (1-NuUtilities::OscillationWeight(energy, dm2, sn2))*
              (fSpectrum->GetBinContent(bin-1) +
               lowGrad*(energy-fSpectrum->GetBinCenter(bin-1)));
            denominator +=
              fSpectrum->GetBinContent(bin-1) +
              lowGrad*(energy-fSpectrum->GetBinCenter(bin-1));
          }
          else {
            oscProb += (1-NuUtilities::OscillationWeight(energy, dm2, sn2))*
              (fSpectrum->GetBinContent(bin) +
               highGrad*(energy-fSpectrum->GetBinCenter(bin)));
            denominator +=
              fSpectrum->GetBinContent(bin) +
              highGrad*(energy-fSpectrum->GetBinCenter(bin));
          }
        }
        if (denominator>0){
          oscProb /= denominator;
        }
        else{
          oscProb = 0;
        }
      }
      
    }
    
    fSpectrum->
      SetBinContent(bin,fSpectrum->GetBinContent(bin)*oscProb);
    fSpectrum->
      SetBinError(bin,fSpectrum->GetBinError(bin)*oscProb);
    
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::InverseOscillateSimpleAverage(const Double_t dm2,
                                                     const Double_t sn2)
{
  for (Int_t bin=1; bin<=fSpectrum->GetNbinsX(); bin++){
    Double_t oscProb=0.0;
    Double_t increment=0.001;
    Int_t count=0;
    for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
         energy<fSpectrum->GetBinLowEdge(bin+1);
         energy+=increment){
      oscProb += 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
      ++count;
    }
    oscProb /= count;
    fSpectrum->SetBinContent(bin,
                             oscProb*fSpectrum->GetBinContent(bin));
    fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::InverseOscillateLinearInterp(const Double_t dm2,
                                                    const Double_t sn2)
{
  for (Int_t bin=1; bin<=fSpectrum->GetNbinsX(); bin++){
    Double_t oscProb=0.0;
    Double_t increment=0.001;
    if (1==bin || fSpectrum->GetNbinsX()==bin){
      Int_t count=0;
      for (Double_t energy = fSpectrum->GetBinLowEdge(bin)+increment/2.0;
           energy<fSpectrum->GetBinLowEdge(bin+1);
           energy+=increment){
        oscProb += 1-NuUtilities::OscillationWeight(energy, dm2, sn2);
        ++count;
      }
      oscProb /= count;
      //fSpectrum->SetBinContent(bin,
      //        oscProb*fSpectrum->GetBinContent(bin));
      //fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
    }
    else{
      Double_t denominator=0.0;
      Double_t lowGrad =
        (fSpectrum->GetBinContent(bin)-fSpectrum->GetBinContent(bin-1))/
        (fSpectrum->GetBinWidth(bin)/2.0 + 
         fSpectrum->GetBinWidth(bin-1)/2.0);
      Double_t highGrad =
        (fSpectrum->GetBinContent(bin+1)-fSpectrum->GetBinContent(bin))/
        (fSpectrum->GetBinWidth(bin+1)/2.0 +
         fSpectrum->GetBinWidth(bin)/2.0);
      for (Double_t energy=fSpectrum->GetBinLowEdge(bin) + increment/2.0;
           energy<fSpectrum->GetBinLowEdge(bin+1);
           energy+=increment){
        if (energy<fSpectrum->GetBinCenter(bin)){
          oscProb += (1-NuUtilities::OscillationWeight(energy, dm2, sn2))*
            (fSpectrum->GetBinContent(bin-1) +
             lowGrad*(energy-fSpectrum->GetBinCenter(bin-1)));
          denominator +=
            fSpectrum->GetBinContent(bin-1) +
            lowGrad*(energy-fSpectrum->GetBinCenter(bin-1));
        }
        else {
          oscProb += (1-NuUtilities::OscillationWeight(energy, dm2, sn2))*
            (fSpectrum->GetBinContent(bin) +
             highGrad*(energy-fSpectrum->GetBinCenter(bin)));
          denominator +=
            fSpectrum->GetBinContent(bin) +
            highGrad*(energy-fSpectrum->GetBinCenter(bin));
        }
      }
      if (denominator>0){
        oscProb /= denominator;
      }
      else{
        oscProb = 0;
      }
    }
  fSpectrum->SetBinContent(bin,oscProb*fSpectrum->GetBinContent(bin));
  fSpectrum->SetBinError(bin, oscProb*fSpectrum->GetBinError(bin));
  }
}

/*
//____________________________________________________________________72
void NuMatrixSpectrum::DecaySimpleAverage(const Double_t alpha,
                                          const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    const Double_t energy = fSpectrum->GetBinCenter(i);
    const Double_t decayProb = NuUtilities::DecayWeight(energy, alpha, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*decayProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*decayProb);
  }
}
*/

//____________________________________________________________________72
void NuMatrixSpectrum::DecayCC(const Double_t alpha,
                             const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    const Double_t energy = fSpectrum->GetBinCenter(i);
    const Double_t decayProb = NuUtilities::DecayWeightCC(energy, alpha, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*decayProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*decayProb);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::DecayNC(const Double_t alpha,
                               const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    const Double_t energy = fSpectrum->GetBinCenter(i);
    const Double_t decayProb = NuUtilities::DecayWeightNC(energy, alpha, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*decayProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*decayProb);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::Decohere(const Double_t mu2,
                                const Double_t sn2)
{
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    const Double_t energy = fSpectrum->GetBinCenter(i);
    const Double_t decoherenceProb = NuUtilities::DecoherenceWeight(energy, mu2, sn2);
    fSpectrum->
      SetBinContent(i,fSpectrum->GetBinContent(i)*decoherenceProb);
    fSpectrum->
      SetBinError(i,fSpectrum->GetBinError(i)*decoherenceProb);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::Add(const TH1D* correction)
{
  // If we don't have a spectrum, copy it
  if (!fSpectrum) {
    fSpectrum = new TH1D(*correction);
  } else {
    // Just add it
    fSpectrum->Add(correction);
  }
}

//____________________________________________________________________72
void NuMatrixSpectrum::Add(const NuMatrixSpectrum& correction, bool addPot)
{
  this->Add(correction.fSpectrum);
  // fSpectrum->Add(correction.fSpectrum);
  if (addPot) fPoT += correction.PoT();
}

//____________________________________________________________________72
void NuMatrixSpectrum::Subtract(const TH1D* correction)
{
  fSpectrum->Add(correction, -1.0);
}

//____________________________________________________________________72
void NuMatrixSpectrum::Subtract(const NuMatrixSpectrum& correction)
{
  fSpectrum->Add(correction.fSpectrum, -1.0);
}


//____________________________________________________________________72
void NuMatrixSpectrum::RebinToArray(std::vector<Double_t> &bins)
{
  TH1D* tmp = NuUtilities::RebinHistogram(fSpectrum, bins);  
  ResetSpectrum(*tmp);
  delete tmp;
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinToArray(int nbins, const double binedges[999])
{
  TH1D* tmp = NuUtilities::RebinHistogram(fSpectrum, nbins, binedges);  
  ResetSpectrum(*tmp);
  delete tmp;
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinToScheme(int scheme)
{
  RebinToScheme(static_cast<NuBinningScheme::NuBinningScheme_t>(scheme));
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinToScheme(NuBinningScheme::NuBinningScheme_t scheme)
{
  TH1D* tmp = NuUtilities::RebinHistogram(fSpectrum, scheme);  
  ResetSpectrum(*tmp);
  delete tmp;  
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinToGeV(Double_t binwidth)
{
  double edges[] = {0, 200};
  double steps[] = {binwidth};
  vector<Double_t> bins = NuUtilities::GenerateBins(1, edges, steps);
  this->RebinToArray(bins);
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinForFit()
{
  double edges[] = {0,   10,  20, 30, 50, 200};
  double steps[] = {0.5, 1.0, 10, 20, 150};
  vector<Double_t> bins = NuUtilities::GenerateBins(5, edges, steps);
  this->RebinToArray(bins);
}


//____________________________________________________________________72
void NuMatrixSpectrum::RebinForPlotAsFit()
{
  double edges[] = {0,   10,  20, 30, 50, 200};
  double steps[] = {0.5, 1.0, 10, 20, 150};
  vector<Double_t> bins = NuUtilities::GenerateBins(5, edges, steps);
  this->RebinToArray(bins);
  this->BinWidthNormalize();
  this->CompressTopBins();
}

//____________________________________________________________________72
void NuMatrixSpectrum::RebinForPlots()
{
  double edges[] = {0,   10,  20, 30, 50, 200};
  double steps[] = {1.0, 2.0, 10, 20, 150};
  vector<Double_t> bins = NuUtilities::GenerateBins(5, edges, steps);
  this->RemoveErrors();
  this->GenerateErrors();
  this->RebinToArray(bins);
  this->BinWidthNormalize();
  this->CompressTopBins();  
}  

//____________________________________________________________________72
NuMatrixSpectrum * NuMatrixSpectrum::Fluctuate()
{
  NuMatrixSpectrum *fluctSpect = new NuMatrixSpectrum(*this);
  fluctSpect->FluctuateMe();
  return fluctSpect;
}

//____________________________________________________________________72
void NuMatrixSpectrum::FluctuateMe()
{
  static NuFluctuator *fl = 0;
  if (!fl) {
    fl = new NuFluctuator();
  }
  
  TH1D tmp = fl->Fluctuate(this->Spectrum());
  ResetSpectrum(tmp);
}

---