NuUtilities.cxx

Go to the documentation of this file.

// Coded by Justin Evans from February 2008 onwards
//
// Contact: justin.evans@ucl.ac.uk

#include <cmath>
#include <glob.h>
#include <sstream>

#include <sys/stat.h>


#include "Rtypes.h"
#include "TFile.h"
#include "TMath.h"

#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "Conventions/ReleaseType.h"
#include "Conventions/SimFlag.h"
#include "MessageService/MsgService.h"
#include "NuMuBar/NuConfig.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuUtilities.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuLibrary.h"

CVSID("$Id: NuUtilities.cxx,v 1.35 2010/02/01 14:41:15 bckhouse Exp $");

//____________________________________________________________________72
NuUtilities::NuUtilities()
{
}

//____________________________________________________________________72
NuUtilities::~NuUtilities()
{
}

//____________________________________________________________________72
void NuUtilities::CalcFiducialMass()
{


  const Float_t avPlSteelDensityNDData=7.847*(Munits::g/Munits::cm3);
  const Float_t avPlSteelThicknessNDData=2.563*(Munits::cm);
  const Float_t avPlScintDensityNDData=1.0457*(Munits::g/Munits::cm3);
  const Float_t avPlScintThicknessNDData=1.02*(Munits::cm);
  const Float_t avPlAlumDensityNDData=2.7*(Munits::g/Munits::cm3);
  const Float_t avPlAlumThicknessNDData=0.1*(Munits::cm);

  const Float_t avPlSteelDensityFDData=avPlSteelDensityNDData;
  const Float_t avPlSteelThicknessFDData=2.559*(Munits::cm);
  const Float_t avPlScintDensityFDData=avPlScintDensityNDData;
  const Float_t avPlScintThicknessFDData=avPlScintThicknessNDData;
  const Float_t avPlAlumDensityFDData=avPlAlumDensityNDData;
  const Float_t avPlAlumThicknessFDData=avPlAlumThicknessNDData;

  //MC
  const Float_t avPlSteelDensityNDMC=avPlSteelDensityNDData;
  const Float_t avPlSteelThicknessNDMC=2.54*(Munits::cm);
  const Float_t avPlScintDensityNDMC=avPlScintDensityNDData;
  const Float_t avPlScintThicknessNDMC=1.0*(Munits::cm);
  const Float_t avPlAlumDensityNDMC=avPlAlumDensityNDData;
  const Float_t avPlAlumThicknessNDMC=avPlAlumThicknessNDData;

  const Float_t avPlSteelDensityFDMC=avPlSteelDensityNDMC;
  const Float_t avPlSteelThicknessFDMC=avPlSteelThicknessNDMC;
  const Float_t avPlScintDensityFDMC=avPlScintDensityNDMC;
  const Float_t avPlScintThicknessFDMC=avPlScintThicknessNDMC;
  const Float_t avPlAlumDensityFDMC=avPlAlumDensityNDMC;
  const Float_t avPlAlumThicknessFDMC=avPlAlumThicknessNDMC;

  cout<<"************** Thickness and Density ******************"<<endl;
  cout<<"Thickness ND Data:"
      <<" steel="<<avPlSteelThicknessNDData
      <<", scint="<<avPlScintThicknessNDData
      <<", alum="<<avPlAlumThicknessNDData<<endl;
  cout<<"Density ND Data:"
      <<" steel="<<avPlSteelDensityNDData
      <<", scint="<<avPlScintDensityNDData
      <<", alum="<<avPlAlumDensityNDData<<endl;
  cout<<endl;
    cout<<"Thickness ND MC:"
      <<" steel="<<avPlSteelThicknessNDMC
      <<", scint="<<avPlScintThicknessNDMC
      <<", alum="<<avPlAlumThicknessNDMC<<endl;
  cout<<"Density ND MC:"
      <<" steel="<<avPlSteelDensityNDMC
      <<", scint="<<avPlScintDensityNDMC
      <<", alum="<<avPlAlumDensityNDMC<<endl;

  cout<<endl<<endl;
  cout<<"Thickness FD Data:"
      <<" steel="<<avPlSteelThicknessFDData
      <<", scint="<<avPlScintThicknessFDData
      <<", alum="<<avPlAlumThicknessFDData<<endl;
  cout<<"Density FD Data:"
      <<" steel="<<avPlSteelDensityFDData
  <<", scint="<<avPlScintDensityFDData
      <<", alum="<<avPlAlumDensityFDData<<endl;
  cout<<endl;
    cout<<"Thickness FD MC:"
      <<" steel="<<avPlSteelThicknessFDMC
      <<", scint="<<avPlScintThicknessFDMC
      <<", alum="<<avPlAlumThicknessFDMC<<endl;
  cout<<"Density FD MC:"
      <<" steel="<<avPlSteelDensityFDMC
      <<", scint="<<avPlScintDensityFDMC
      <<", alum="<<avPlAlumDensityFDMC<<endl;

  const Float_t avPlSteelStPwNDData=
    avPlSteelThicknessNDData*avPlSteelDensityNDData;
  const Float_t avPlScintStPwNDData=
    avPlScintThicknessNDData*avPlScintDensityNDData;
  const Float_t avPlAlumStPwNDData=
    avPlAlumThicknessNDData*avPlAlumDensityNDData;

  const Float_t avPlSteelStPwFDData=
    avPlSteelThicknessFDData*avPlSteelDensityFDData;
  const Float_t avPlScintStPwFDData=
    avPlScintThicknessFDData*avPlScintDensityFDData;
  const Float_t avPlAlumStPwFDData=
    avPlAlumThicknessFDData*avPlAlumDensityFDData;

  const Float_t avPlSteelStPwNDMC=
    avPlSteelThicknessNDMC*avPlSteelDensityNDMC;
  const Float_t avPlScintStPwNDMC=
    avPlScintThicknessNDMC*avPlScintDensityNDMC;
  const Float_t avPlAlumStPwNDMC=
    avPlAlumThicknessNDMC*avPlAlumDensityNDMC;

  const Float_t avPlSteelStPwFDMC=
    avPlSteelThicknessFDMC*avPlSteelDensityFDMC;
  const Float_t avPlScintStPwFDMC=
    avPlScintThicknessFDMC*avPlScintDensityFDMC;
  const Float_t avPlAlumStPwFDMC=
    avPlAlumThicknessFDMC*avPlAlumDensityFDMC;

  cout<<endl<<endl;
  cout<<"**************** Stopping Power ********************"<<endl;
  cout<<"StPw ND Data:"
      <<" steel="<<avPlSteelStPwNDData
      <<", scint="<<avPlScintStPwNDData
      <<", alum="<<avPlAlumStPwNDData<<endl;
  cout<<endl;
  cout<<"StPw ND MC:"
      <<" steel="<<avPlSteelStPwNDMC
      <<", scint="<<avPlScintStPwNDMC
      <<", alum="<<avPlAlumStPwNDMC<<endl;
  cout<<endl<<endl;
  cout<<"StPw FD Data:"
      <<" steel="<<avPlSteelStPwFDData
      <<", scint="<<avPlScintStPwFDData
      <<", alum="<<avPlAlumStPwFDData<<endl;
  cout<<endl;
  cout<<"StPw FD MC:"
      <<" steel="<<avPlSteelStPwFDMC
      <<", scint="<<avPlScintStPwFDMC
      <<", alum="<<avPlAlumStPwFDMC<<endl;


  //NOW CALC AREAS
  Float_t areaPlND=3.142*pow(0.8*Munits::m,2);
  Float_t areaPlFD=3.142*(14*Munits::m2-pow(0.4*Munits::m,2));
  cout<<endl<<endl;
  cout<<"Area ND = "<<areaPlND<<endl;
  cout<<"Area FD = "<<areaPlFD<<endl;

  //NOW CALC MASS
  const Float_t avPlSteelMassNDData=avPlSteelStPwNDData*areaPlND;
  const Float_t avPlScintMassNDData=avPlScintStPwNDData*areaPlND;
  const Float_t avPlAlumMassNDData=avPlAlumStPwNDData*areaPlND;

  const Float_t avPlSteelMassNDMC=avPlSteelStPwNDMC*areaPlND;
  const Float_t avPlScintMassNDMC=avPlScintStPwNDMC*areaPlND;
  const Float_t avPlAlumMassNDMC=avPlAlumStPwNDMC*areaPlND;

  const Float_t avPlSteelMassFDData=avPlSteelStPwFDData*areaPlFD;
  const Float_t avPlScintMassFDData=avPlScintStPwFDData*areaPlFD;
  const Float_t avPlAlumMassFDData=avPlAlumStPwFDData*areaPlFD;

  const Float_t avPlSteelMassFDMC=avPlSteelStPwFDMC*areaPlFD;
  const Float_t avPlScintMassFDMC=avPlScintStPwFDMC*areaPlFD;
  const Float_t avPlAlumMassFDMC=avPlAlumStPwFDMC*areaPlFD;

  cout<<endl<<endl;
  cout<<"**************** MASS ********************"<<endl;
  cout<<"Mass ND Data:"
      <<" steel="<<avPlSteelMassNDData
      <<", scint="<<avPlScintMassNDData
      <<", alum="<<avPlAlumMassNDData<<endl;
  cout<<endl;
  cout<<"Mass ND MC:"
      <<" steel="<<avPlSteelMassNDMC
      <<", scint="<<avPlScintMassNDMC
      <<", alum="<<avPlAlumMassNDMC<<endl;
  cout<<endl<<endl;
  cout<<"Mass FD Data:"
      <<" steel="<<avPlSteelMassFDData
      <<", scint="<<avPlScintMassFDData
      <<", alum="<<avPlAlumMassFDData<<endl;
  cout<<endl;
  cout<<"Mass FD MC:"
      <<" steel="<<avPlSteelMassFDMC
      <<", scint="<<avPlScintMassFDMC
      <<", alum="<<avPlAlumMassFDMC<<endl;



  const Float_t avPlTotalMassNDData=
    avPlSteelMassNDData+avPlScintMassNDData+avPlAlumMassNDData;
  const Float_t avPlTotalMassNDMC=
    avPlSteelMassNDMC+avPlScintMassNDMC+avPlAlumMassNDMC;
  const Float_t avPlTotalMassFDData=
    avPlSteelMassFDData+avPlScintMassFDData+avPlAlumMassFDData;
  const Float_t avPlTotalMassFDMC=
    avPlSteelMassFDMC+avPlScintMassFDMC+avPlAlumMassFDMC;

  const Float_t massSumND=avPlTotalMassNDData+avPlTotalMassNDMC;
  const Float_t massSumFD=avPlTotalMassFDData+avPlTotalMassFDMC;
  const Float_t massDiffND=avPlTotalMassNDData-avPlTotalMassNDMC;
  const Float_t massDiffFD=avPlTotalMassFDData-avPlTotalMassFDMC;

  const Float_t massRelDiffND=massDiffND/(0.5*massSumND);
  const Float_t massRelDiffFD=massDiffFD/(0.5*massSumFD);

  cout<<endl<<endl;
  cout<<"Mass ND Data Total: "<<avPlTotalMassNDData<<endl;
  cout<<"Mass ND MC Total: "<<avPlTotalMassNDMC<<endl;
  cout<<"Mass ND Data-MC="<<massDiffND<<endl;
  cout<<"Mass ND Data-MC/(0.5*Data+MC)="<<massRelDiffND<<endl;
  cout<<endl;
  cout<<"Mass FD Data Total: "<<avPlTotalMassFDData<<endl;
  cout<<"Mass FD MC Total: "<<avPlTotalMassFDMC<<endl;
  cout<<"Mass FD Data-MC="<<massDiffFD<<endl;
  cout<<"Mass FD Data-MC/(0.5*Data+MC)="<<massRelDiffFD<<endl;
  cout<<endl;

  //Considering events from STEEL planes
  //ND: 14->68 inclusive = 55 planes of steel
  cout<<endl<<endl;
  cout<<"55 planes in ND Data="
      <<55*avPlTotalMassNDData/(Munits::kilotonne)<<" kilotonnes"<<endl;
  cout<<"55 planes in ND MC  ="
      <<55*avPlTotalMassNDMC/(Munits::kilotonne)<<" kilotonnes"<<endl;

  //FD: 4->239 inclusive + 253->464 inclusive = 236 + 212 = 448 planes
  cout<<endl<<endl;
  cout<<"448 planes in FD Data="
      <<448*avPlTotalMassFDData/(Munits::kilotonne)<<" kilotonnes"
      <<endl;
  cout<<"448 planes in FD MC  ="
      <<448*avPlTotalMassFDMC/(Munits::kilotonne)<<" kilotonnes"
      <<endl;


//Table 5 of minos-doc-3134v2 has updated numbers
//Both detectors have the same densities in MC
//steel = 7.847
//scint = 1.0457, thickness=1.0cm
//The pitches of ND and FD are slightly different

//Peter Litchfield's email of 20/06/07 22:27 GMT
//Real detectors
//Plane pitch  FD (average) 5.946cm,  ND 5.94cm (3134)
//Steel density  Both Detectors 7.847gm/cc (1431)
//Steel thickness FD 2.559cm, ND 2.563cm (1529)
//Scintillator 1.02cm density 1.0457gm/cc (scint+coating, 3134)
//Aluminium    0.1cm  density 2.7gm/cc (PDG book)
//FD average density = 3.602gm/cc
//ND average density = 3.611gm/cc

//In the Birch MC (3134)
//Plane pitch  FD (average) 5.946cm,  ND 5.94cm
//Steel density 7.87gm/cc
//Steel thickness 2.54cm
//Scintillator 1.02cm density 1.032gm/cc
//Aluminium 0.1cm air, negligible density
//FD avearge density = 3.539gm/cc
//ND average density = 3.542gm/cc
}

//____________________________________________________________________72
TString NuUtilities::PrintRelease(int releaseType)
{
  char hexnum[20];
  sprintf(hexnum, "%x", releaseType);
  TString ret = ReleaseType::AsString(releaseType);
  ret += " (0x";
  ret += hexnum;
  ret += ", ";
  ret += releaseType;
  ret += ")";
  return ret;
}


//____________________________________________________________________72
const BeamType::BeamType_t NuUtilities::ReturnConventionsBeamType
(const NuConfig& config)
{
  if (Detector::kFar==config.detector && SimFlag::kData==config.simFlag){
    MAXMSG("NuUtilities",Msg::kInfo,1)
      << "Detected FD data, so not combining intensity with beamType"
      << endl;
    BeamType::BeamType_t originalBeamType =
      static_cast<BeamType::BeamType_t>(config.beamType);
    return originalBeamType;
  }
  else{
    return NuUtilities::ConvertToConventionsBeamType(config.beamType,
                                                     config.intensity);
  }
}

//____________________________________________________________________72
const BeamType::BeamType_t NuUtilities::ConvertToConventionsBeamType
(const Int_t beamType, const Int_t intensity)
{
  BeamType::BeamType_t originalBeamType =
    static_cast<BeamType::BeamType_t>(beamType);
  if (0==intensity){
    return originalBeamType;
  }
  if (BeamType::kL010z185i == originalBeamType){
    if (124==intensity){
      return BeamType::kL010z185i_i124;
    }
    else if (191==intensity){
      return BeamType::kL010z185i_i191;
    }
    else if (213==intensity){
      return BeamType::kL010z185i_i213;
    }
    else if (224==intensity){
      return BeamType::kL010z185i_i224;
    }
    else if (232==intensity){
      return BeamType::kL010z185i_i232;
    }
    else if (243==intensity){
      return BeamType::kL010z185i_i243;
    }
    else if (257==intensity){
      return BeamType::kL010z185i_i257;
    }
    else if (282==intensity){
      return BeamType::kL010z185i_i282;
    }
    else if (303==intensity){
      return BeamType::kL010z185i_i303;
    }
    else if (324==intensity){
      return BeamType::kL010z185i_i324;
    }
    else {
      MSG("NuUtilities",Msg::kWarning)
        << "Received an L010z185i beam type with unrecognized intensity. Intensity = "
        << intensity
        << endl;
      return BeamType::kUnknown;
    }
  }

  if (BeamType::kL010z000i == originalBeamType){
    if (209 == intensity){
      return BeamType::kL010z000i_i209;
    }
    else if (225 == intensity){
      return BeamType::kL010z000i_i225;
    }
    else if (232 == intensity){
      return BeamType::kL010z000i_i232;
    }
    else if (259 == intensity){
      return BeamType::kL010z000i_i259;
    }
    else if (300 == intensity){
      return BeamType::kL010z000i_i300;
    }
    else if (317 == intensity){
      return BeamType::kL010z000i_i317;
    }
    else if (326 == intensity){
      return BeamType::kL010z000i_i326;
    }
    else if (380 == intensity){
      return BeamType::kL010z000i_i380;
    }
    else {
      MSG("NuUtilities",Msg::kWarning)
        << "Received an L010z000i beam type with unrecognized intensity."
        << endl;
      return BeamType::kUnknown;
    }
  }

  if (BeamType::kL250z200i == originalBeamType){
    if (100 == intensity){
      return BeamType::kL250z200i_i100;
    }
    else if (114 == intensity){
      return BeamType::kL250z200i_i114;
    }
    else if (130 == intensity){
      return BeamType::kL250z200i_i130;
    }
    else if (152 == intensity){
      return BeamType::kL250z200i_i152;
    }
    else if (165 == intensity){
      return BeamType::kL250z200i_i165;
    }
    else if (194 == intensity){
      return BeamType::kL250z200i_i194;
    }
    else if (232 == intensity){
      return BeamType::kL250z200i_i232;
    }
    else {
      MSG("NuUtilities",Msg::kWarning)
        << "Received an L250z200i beam type with unrecognized intensity."
        << endl;
      return BeamType::kUnknown;
    }
  }

  MSG("NuUtilities",Msg::kWarning)
    << "Received a non-zero intensity flag, along with a beam "
    << "type not configured for intensity weighting."
    << endl;
  return BeamType::kUnknown;

}

//____________________________________________________________________72
const vector<Double_t> NuUtilities::GenerateBins(int ngroups, double edges[], double steps[])
{
    vector<Double_t> bins;
    for (int g = 0; g < ngroups; g++) {
        Double_t lowEdge = edges[g];
        Double_t highEdge = edges[g+1];
        Double_t binWidth = steps[g];
        for (Double_t binEdge = lowEdge;
             binEdge < highEdge - (binWidth/2.0);
             binEdge += binWidth){
            bins.push_back(binEdge);
        }
    }
    bins.push_back(edges[ngroups]);
    return bins;
}


//____________________________________________________________________72
const vector<Double_t> NuUtilities::RecoBins
(const NuBinningScheme::NuBinningScheme_t binningScheme)
{
  if (NuBinningScheme::kUnknown == binningScheme ||
      NuBinningScheme::k0_5GeVTo200 == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "Binning reco energy in 0.5 GeV up to 200 GeV" << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
  else if (NuBinningScheme::kCC0325Std == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "Reco energy: one bin from 0.0--0.5, "
      << "then 0.25 GeV bins up to 200.0 GeV"
      << endl;
    vector<Double_t> bins;
    bins.push_back(0.0);
    Double_t lowEdge = 0.5;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.25;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
  else if (NuBinningScheme::kNuMuBar0325Std == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "Reco energy: 0.5 GeV up to 30 GeV; 2 GeV up to 50 GeV; "
      << "1 bin up to 200 GeV"
      << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 30.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge - (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge - (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kNuMuBar0325Std2 == binningScheme)
  {
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "Reco Energy: One 0.5 GeV, then 0.25 GeV up to 20 GeV; "
      << "1 Gev to 30 GeV,\n 2 GeV from 30 to 50 GeV; One 150 GeV bin to 200 GeV"
      << endl;
    vector<Double_t> bins;
    bins.push_back(0.0);
    Double_t lowEdge = 0.5;
    Double_t highEdge = 20.0;
    Double_t binWidth = 0.25;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 20.0;
    highEdge = 30.0;
    binWidth = 1.0;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kDisplayFD == binningScheme) {
    MAXMSG("NuUtilities",Msg::kInfo,5)
    << "Reco Energy: 4 GeV Bins up to 20 GeV, then 20, 30, 50, 200"
    << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 20.0;
    Double_t binWidth = 4;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(20.0);
    bins.push_back(30.0);
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kDisplayND == binningScheme) {
    MAXMSG("NuUtilities",Msg::kInfo,5)
    << "Reco Energy: 1 Gev to 30 GeV, 2 GeV from 30 to 50 GeV; One 150 GeV bin to 200 GeV"
    << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 20.0;
    Double_t binWidth = 1.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 20.0;
    highEdge = 30.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 4.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else {
    MAXMSG("NuUtilities",Msg::kInfo,100)
    << "!!Invalid binning scheme (" << binningScheme << ") requested.!!" << endl
    << "Binning reco energy in 0.5 GeV up to 200 GeV" << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
}

//____________________________________________________________________72
const vector<Double_t> NuUtilities::TrueBins
(const NuBinningScheme::NuBinningScheme_t binningScheme)
{
  if (NuBinningScheme::kUnknown == binningScheme ||
      NuBinningScheme::k0_5GeVTo200 == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "Binning true energy in 0.5 GeV up to 200 GeV" << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
  else if (NuBinningScheme::kCC0325Std == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "True energy: one bin from 0.0--0.5, "
      << "then 0.25 GeV bins up to 200.0 GeV"
      << endl;
    vector<Double_t> bins;
    bins.push_back(0.0);
    Double_t lowEdge = 0.5;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.25;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
  else if (NuBinningScheme::kNuMuBar0325Std == binningScheme){
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "True energy: 0.5 GeV up to 30 GeV; 2 GeV up to 50 GeV; "
      << "1 bin up to 200 GeV"
      << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 30.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge - (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge - (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kNuMuBar0325Std2 == binningScheme) {
    MAXMSG("NuUtilities",Msg::kInfo,5)
      << "True Energy: One 0.5 GeV, then 0.25 GeV up to 20 GeV; "
      << "1 Gev to 30 GeV,\n 2 GeV from 30 to 50 GeV; One 150 GeV bin to 200 GeV"
      << endl;
    vector<Double_t> bins;
    bins.push_back(0.0);
    Double_t lowEdge = 0.5;
    Double_t highEdge = 20.0;
    Double_t binWidth = 0.25;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 20.0;
    highEdge = 30.0;
    binWidth = 1.0;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
          binEdge < highEdge -(binWidth/2.0);
          binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kDisplayFD == binningScheme) {
    MAXMSG("NuUtilities",Msg::kInfo,5)
    << "True Energy: 4 GeV Bins up to 20 GeV, then 20, 30, 50, 200"
    << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 20.0;
    Double_t binWidth = 4;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(20.0);
    bins.push_back(30.0);
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else if (NuBinningScheme::kDisplayND == binningScheme) {
    MAXMSG("NuUtilities",Msg::kInfo,5)
    << "Reco Energy: 1 Gev to 30 GeV, 2 GeV from 30 to 50 GeV; One 150 GeV bin to 200 GeV"
    << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 20.0;
    Double_t binWidth = 1.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 20.0;
    highEdge = 30.0;
    binWidth = 2.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    lowEdge = 30.0;
    highEdge = 50.0;
    binWidth = 4.0;
    for (Double_t binEdge = lowEdge;
         binEdge < highEdge -(binWidth/2.0);
         binEdge += binWidth) {
      bins.push_back(binEdge);
    }
    bins.push_back(50.0);
    bins.push_back(200.0);
    return bins;
  }
  else {
    MAXMSG("NuUtilities",Msg::kInfo,100)
    << "!!Invalid binning scheme requested.!!" << endl
    << "Binning true energy in 0.5 GeV up to 200 GeV" << endl;
    vector<Double_t> bins;
    Double_t lowEdge = 0.0;
    Double_t highEdge = 200.0;
    Double_t binWidth = 0.5;
    for (Double_t binEdge = lowEdge;
         binEdge <= highEdge + (binWidth/2.0);
         binEdge += binWidth){
      bins.push_back(binEdge);
    }
    return bins;
  }
}

//____________________________________________________________________72
void NuUtilities::WriteHistosToFile(const std::vector<TH1D>& vHistos,
                                    TFile& file)
{
  file.cd();
  for (vector<TH1D>::const_iterator it = vHistos.begin();
       it != vHistos.end();
       ++it){
    (*it).Write();
    cout << "Written " << (*it).GetName() << endl;
  }
}

//____________________________________________________________________72
void NuUtilities::WriteHistosToFile(const std::vector<NuMatrixSpectrum>& vHistos,
                                    TFile& file)
{
  file.cd();
  for (vector<NuMatrixSpectrum>::const_iterator it = vHistos.begin();
       it != vHistos.end();
       ++it){
    (*it).Spectrum()->Write();
    cout << "Written " << (*it).GetName() << endl;
  }
}


//____________________________________________________________________72
TH1D* NuUtilities::RebinHistogram(const TH1D* old, int nbins, const double binedges[999])
{
  static int count = 0;
  TString scount = "";
  scount += count++;

  TH1D *newh = new TH1D(old->GetName()+TString("_rebin")+scount, old->GetTitle(), nbins, binedges);
  bool ouwarn = true;

  // Now, loop over the old histogram and fill the new one
  for(Int_t bin = 1; bin <= old->GetNbinsX(); bin++) {

    double oldCent = old->GetBinLowEdge(bin);
    int newbin = newh->GetXaxis()->FindFixBin(oldCent);

    // Make sure that bins do not cross boundaries.
    if (newbin > newh->GetNbinsX() || newbin == 0) {
      if (ouwarn) {
        cout << "Old bin: " << old->GetBinLowEdge(bin) << " - " << old->GetBinLowEdge(bin+1)
        << " and maybe others have become under/overflow" << endl;
        ouwarn = false;
      }
    }
    else if ( old->GetBinLowEdge(bin) < newh->GetBinLowEdge(newbin) ||
        old->GetBinLowEdge(bin+1) > newh->GetBinLowEdge(newbin+1)) {
      cout << "Bin edge conflict in bin no's " << bin << ", " << newbin << ".  No rebinning performed."
           << "  old bin: " << old->GetBinLowEdge(bin) << " - " << old->GetBinLowEdge(bin+1)
           << ", new bin: " << newh->GetBinLowEdge(newbin) << " - " << newh->GetBinLowEdge(newbin+1)
           << endl;
      TH1D *old_ret = new TH1D(*old);
      return old_ret;
    }

    double oldCon = old->GetBinContent(bin);
    double oldErr = old->GetBinError(bin);
    double newCon = newh->GetBinContent(newbin);
    double newErr = newh->GetBinError(newbin);

    newCon += oldCon;
    newErr = sqrt(newErr*newErr + oldErr*oldErr);
    newh->SetBinContent(newbin, newCon);
    newh->SetBinError(newbin, newErr);
  }

  // Validate
  double oldint = old->Integral();
  double newint = newh->Integral() + newh->GetBinContent(0)
                 + newh->GetBinContent(newh->GetNbinsX()+1);
  if (abs(oldint - newint) > 1e-5) {
    cout << "Rebinning report: " << oldint << " entries have become " << newint << " events" << endl;
  }
  //newh->Sumw2();
  return newh;
}


//____________________________________________________________________72
TH1D* NuUtilities::RebinHistogram(const TH1D* old, std::vector<Double_t> &bins)
{
  return RebinHistogram(old, bins.size()-1, &(bins[0]));
}

//____________________________________________________________________72
TH1D* NuUtilities::RebinHistogram(const TH1D* old, NuBinningScheme::NuBinningScheme_t scheme)
{
  vector<Double_t> bins = RecoBins(scheme);
  return RebinHistogram(old, bins);
}


template<class T> inline T SQR(T x){return x*x;}


//____________________________________________________________________72
Double_t NuUtilities::OscillationWeight(const Double_t energy,
                                        const Double_t dm2,
                                        const Double_t sn2)
{
  return 1-sn2*SQR(sin(1.27*dm2*(FDDistance()/Munits::km)/energy));
}


//____________________________________________________________________72
Double_t NuUtilities::DecayWeightCC(const Double_t energy,
                                  const Double_t alpha,
                                  const Double_t sn2)
{
  return SQR(sn2+(1-sn2)*exp(-(alpha*FDDistance()/Munits::km)/(2*energy)));
}


//____________________________________________________________________72
Double_t NuUtilities::DecayWeightNC(const Double_t energy,
                                    const Double_t alpha,
                                    const Double_t sn2)
{
  return (1-sn2)+sn2*exp(-(alpha*FDDistance()/Munits::km)/energy);
}


//____________________________________________________________________72
Double_t NuUtilities::DecoherenceWeight(const Double_t energy,
                                        const Double_t mu2,
                                        const Double_t sn2)
{
  return 1-(sn2/2)*(1-exp(-(mu2*FDDistance()/Munits::km)/(2*energy)));
}

//____________________________________________________________________72
void NuUtilities::FixDogwoodQP(NuEvent &nu)
{
  // Only need to fix dogwood NuEvent's
  if (!ReleaseType::IsDogwood(nu.recoVersion)) {
    MAXMSG("NuUtilities",Msg::kInfo,5) << "Not correcting q/p for "
    << PrintRelease(nu.recoVersion) << endl;
    return;
  }
  MAXMSG("NuUtilities",Msg::kInfo,5) << "Correcting "
  << PrintRelease(nu.recoVersion) << " q/p from "
  << nu.qp << " to " << nu.qp_rangebiased << endl;

  nu.qp = nu.qp_rangebiased;
  if (nu.sigqp == 0) nu.qp_sigqp = 0;
  else               nu.qp_sigqp = nu.qp_rangebiased / nu.sigqp;

  nu.qp1 = nu.qp_rangebiased1;
  if (nu.sigqp1 == 0) nu.qp_sigqp1 = 0;
  else                nu.qp_sigqp1 = nu.qp_rangebiased1 / nu.sigqp1;

  nu.qp2 = nu.qp_rangebiased2;
  if (nu.sigqp2 == 0) nu.qp_sigqp2 = 0;
  else                nu.qp_sigqp2 = nu.qp_rangebiased2 / nu.sigqp2;

  nu.qp3 = nu.qp_rangebiased3;
  if (nu.sigqp3 == 0) nu.qp_sigqp3 = 0;
  else                nu.qp_sigqp3 = nu.qp_rangebiased3 / nu.sigqp3;

  if (nu.qp_sigqp == 0) nu.sigqp_qp = 0;
  else                  nu.sigqp_qp = 1./nu.qp_sigqp;

  NuLibrary& lib=NuLibrary::Instance();
  lib.reco.GetEvtEnergy(nu, false);

  if (nu.qp_rangebiased > 0) nu.charge = 1;
  else                       nu.charge = -1;
}

//____________________________________________________________________72
std::vector<TString> NuUtilities::
  GetListOfFilesInDir(const TString &wildcardString)
{
  std::vector<TString> fileList;

  glob_t g;
  glob(wildcardString.Data(),
       // GLOB_NOCHECK | // If no match return pattern
       GLOB_TILDE, // Expand ~'s
       0, &g);
  for(unsigned int i = 0; i < g.gl_pathc; ++i)
    fileList.push_back(g.gl_pathv[i]);
  globfree(&g);

  // glob output is sorted by default

  return fileList;
}
//____________________________________________________________________72

Double_t NuUtilities::LogLikelihood(Double_t dnu, Double_t mnu)
{
  if (dnu) {
    // Protection against cases where there is data but no MC?
    // This should not normally be possible, other than highly contrived
    // unphysical oscillation values.
    // Set the value to something lower than POT_Data / POT_MC, which at the
    // moment is usually 1e-3. We can justify this as an error on the MC
    const Double_t MCError = 1e-7;
    if (mnu < MCError) mnu = MCError;
        
    // We have both values. Do the proper lnL!
    return 2*(mnu - dnu + dnu*log(dnu/mnu));
  } else {
    // dnu is zero, mnu may or may not be zero
    return 2*mnu;
  }
}


//____________________________________________________________________72

int NuUtilities::GetDigit(int num, int d) 
{ 
  int pow1 = static_cast<int>(TMath::Power(10, d));
  int pow2 = pow1 / 10;
  return (num % pow1 - num % pow2) / pow2; 
}



//____________________________________________________________________72

void NuUtilities::ProgressBar(Int_t entry, Float_t nEntriesF, Int_t mintime)
{
  Int_t nEntries = (Int_t)nEntriesF;
    
  // Are we streaming to a file? If so, show the old style
  struct stat buf;
  fstat(fileno(stdout), &buf);
  // Check if we are a file, or a pipe (i.e. in case the user tees)
  const bool isFile = buf.st_mode & (S_IFREG | S_IFIFO) ;
  if (isFile) {
    Float_t fract = ceil(nEntries/20.);  
    if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract)
    {
        MSG("NuUtilities",Msg::kInfo) 
          <<"Fraction of loop complete: "<<entry 
          <<"/"<<nEntries<<"  ("
          <<(Int_t)(100.*entry/nEntries)<<"%)"<<endl;
    }
    return;
  }
  
  
  // How wide should we make the progress bar?
  const Int_t width = 70;
  // How long is the string for entries?
  static Int_t countlen = -1;
  // How long is our progress bar?
  static Int_t barlen = -1;
  // The entry number of the next bar entry
  static Int_t nextbar = -1;
  // When did we start?
  static time_t starttime = 0;
  // when do we next update?
  static time_t nextupdate = 0;
  // If we are processing the first entry, reset everything
  if (entry <= 1)
  {
    // Get the new length of the entry string
    countlen = (Int_t)TMath::Ceil(TMath::Log10(nEntries)) + 1;
    nextbar = -1;
    starttime = time(NULL);

    barlen = width - 14 - countlen*2 - 1;

    // Don't update until we get to the minimum time
    nextupdate = starttime + mintime;
  }

  // Check here to see if we should update; otherwise, return
  // Check to see if the bar would update
  // or, alternatively, it is time to refresh.
  if ((time(NULL) < nextupdate) && (entry < nextbar || (nextbar == -1)))return;
  nextupdate = time(NULL) + 10;

  // Because this is used in several places, make it here
  Float_t frac = (Float_t)entry / (Float_t)nEntries;

  // Prepare the progress bar string
  TString bar;
  if (entry <= nEntries)
  {
    // Work out how many characters we are in
    Int_t numeq = TMath::FloorNint(frac*barlen);

    // Work out when the next bar will occur
    nextbar = (Int_t)((Float_t)(numeq+1)/(Float_t)barlen*nEntries);
    //cout << "Next bar at: " << nextbar << "        " << endl;
    bar = TString('=', numeq);
    bar += TString(' ', barlen - numeq);
  } else if (entry > nEntries) {
    // We have gone over. Oh no!
    bar = TString('+', barlen);
  } else if (entry < 0) {
    // Somehow, we are below zero. Handle it nonetheless
    bar = TString('-', barlen);
  }


  // Prepare the ETA
  Float_t elapsed_time = (Float_t)(time(NULL) - starttime);
  Float_t time_left = -60;
  if (frac > 1e-6) {
    time_left = (elapsed_time / frac) - elapsed_time;
  }
  Int_t mins, seconds;
  mins    = (Int_t)TMath::Floor(time_left / 60.0f);
  seconds = (Int_t)TMath::Floor(time_left - (Float_t)(mins*60.0f));
  // TString ETA;
  ostringstream ETA;

  ETA << "ETA ";
  if ((mins < 0 || seconds < 0) || (mins == 0 && seconds == 0)) {
    ETA << "--:--";
  } else {
    ETA << setfill('0') << setw(2) << mins << ":" << setw(2) << seconds;
  }

  cout << " Progress: [" << bar << "] "
       << setw(countlen) << entry << "/" << nEntries
       << " " << ETA.str()
       << '\r'
       << flush;
  // Move to the next line, if this is the final entry!
  if (entry == nEntries) {
    cout << endl;
  }
}

---