NuMatrix1D.cxx

Go to the documentation of this file.

// NuMatrix1D.cxx
// A slightly less abstract, 1D data storage class

#include <cassert>

#include "NuMatrix1D.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuUtilities.h"

#include "MessageService/MsgService.h"
#include "Conventions/Munits.h"

#include "TFile.h"
#include "TH1D.h"
#include "TMath.h"
#include "TGraph.h"

CVSID("$Id: NuMatrix1D.cxx,v 1.3 2010/01/11 09:43:30 nickd Exp $");

ClassImp(NuMatrix1D)


// Constructors

NuMatrix1D::NuMatrix1D() : NuMatrix(), fSpectrum(0)
{
  
}

// The POT-only constructor
NuMatrix1D::NuMatrix1D(Double_t POT) : NuMatrix(POT), fSpectrum(0)
{
  
}

// The spectrum-only constructor
NuMatrix1D::NuMatrix1D(const TH1D& spectrum) : NuMatrix(0), fSpectrum(0)
{
  fSpectrum = new TH1D(spectrum);
  fSpectrum->SetDirectory(0);
}

// The spectrum *and* POT constructor
NuMatrix1D::NuMatrix1D(const TH1D& spectrum, Double_t POT) : NuMatrix(POT), fSpectrum(0)
{
  fSpectrum = new TH1D(spectrum);
  fSpectrum->SetDirectory(0);
}

// The copy-constructor
NuMatrix1D::NuMatrix1D(const NuMatrix1D& original) : NuMatrix(original), fSpectrum(0)
{
  if (original.fSpectrum) {
    fSpectrum = new TH1D(*original.fSpectrum);
    fSpectrum->SetDirectory(0);
  }
}

// The assignment operator
NuMatrix1D & NuMatrix1D::operator=(const NuMatrix1D& orig)
{
  if (this != &orig) {
    // Copy the base class stuff
    this->NuMatrix::operator=(orig);
  
    // Copy the histogram
    if (fSpectrum) delete fSpectrum; fSpectrum = 0;
    if (orig.fSpectrum) {
      fSpectrum = new TH1D(*orig.fSpectrum);
      fSpectrum->SetDirectory(0);
    }
  }  
  return *this;
}

// Copy a NuMatrixSpectrum
NuMatrix1D::NuMatrix1D(const NuMatrixSpectrum& original) : fSpectrum(0)
{
  // Copy the POT from the NuMatrixSpectrum
  ResetPOT(original.PoT());
  
  // Copy the spectrum from a NuMatrixSpectrum!
  fSpectrum = new TH1D(*original.Spectrum());
}

// Assigned from a NuMatrixspectrum
NuMatrix1D &NuMatrix1D::operator=(const NuMatrixSpectrum& orig)
{
  // Copy the POT count
  this->ResetPOT(orig.PoT());
  
  // Copy the spectrum from a NuMatrixSpectrum!
  if (fSpectrum) delete fSpectrum;
  fSpectrum = new TH1D(*orig.Spectrum());
  
  return *this;
}


NuMatrix1D::~NuMatrix1D()
{
  if (fSpectrum) delete fSpectrum; fSpectrum = 0;
}

// Regular member functions ////////////////////////////////////////////

Double_t NuMatrix1D::StatsLikelihood(const NuMatrix* with) const
{
  // Convert the NuMatrix to a NuMatrix1D
  const NuMatrix1D *other = dynamic_cast<const NuMatrix1D*>(with);
  if (other == 0) {
    MSG("NuMatrix1D",Msg::kWarning) << "Attempting to compare a NuMatrix1D to something other than a NuMatrix1D" << endl;
    return -1;
  }
  
  // // Do the statistical comparison
  // MSG("NuMatrix1D",Msg::kError) << "Comparing statistical likelihood - but function not implemented" << endl;
  // return 0;
  return this->StatsLikelihood(other);
}


void NuMatrix1D::ScaleToPOT(Double_t new_POT)
{
  // check we have POT and a spectrum
  if (GetPOT() != 0.0 && fSpectrum) {
    Double_t scale_factor = new_POT / GetPOT();
    fSpectrum->Scale(scale_factor);
  }
  
  // Reset the internal POT
  ResetPOT(new_POT);
}


// Load the 1D Matrix from a histogram, with a given POT
NuMatrix1D::NuMatrix1D( const TString &filename,
                        const TString &histname,
                        Double_t POT)
  : NuMatrix(), fSpectrum(0)
{
  // Attempt to open the named file
  TFile tf(filename);
  if (!tf.IsOpen())
  {
    MSG("NuMatrix1D",Msg::kError) << "Cannot open file " << filename << endl;
    return;
  }
  
  // Now we have opened, grab the histogram
  fSpectrum = dynamic_cast<TH1D*>(tf.Get(histname));
  // Detach the spectrum from the file
  fSpectrum->SetDirectory(0);
  
  // Check this worked
  if (fSpectrum == 0)
  {
    MSG("NuMatrix1D",Msg::kError) << "Cannot read histogram " << histname
      << " from file " << filename << endl;
    return;
  }
  
  // Now we have successfully read the histogram, set the POT
  this->ResetPOT(POT);
  
}


// Load the 1D matrix from a file, and read the POT from another histogram
NuMatrix1D::NuMatrix1D( const TString &filename,
                        const TString &histname,
                        const TString &pothistname)
                      : NuMatrix(), fSpectrum(0)
{
  // Attempt to open the named file
  TFile tf(filename);
  if (!tf.IsOpen())
  {
    MSG("NuMatrix1D",Msg::kError) << "Cannot open file " << filename << endl;
    return;
  }
  
  // Only try to read the histogram if we have been given a name!
  if (histname != "") {
    // Now we have opened, grab the histogram
    TH1D *spect = dynamic_cast<TH1D*>(tf.Get(histname));
    // Check this worked
    if (spect == 0)
    {
      MSG("NuMatrix1D",Msg::kError) << "Cannot read histogram " << histname
        << " from file " << filename << endl;
      return;
    }
    fSpectrum = new TH1D(*spect);
    // Detach the spectrum from the file
    fSpectrum->SetDirectory(0);
  }
  
  // Now we have successfully read the histogram, get the POT histogram
  TH1 *hTotalPot = dynamic_cast<TH1*>(tf.Get(pothistname));
  if (hTotalPot == 0)
  {
    MSG("NuMatrix1D",Msg::kError) << "Cannot read POT histogram " << pothistname
      << " from file " << filename << endl;
    fSpectrum = 0;
    return;
  }
  // Set the POT from this histogram
  this->ResetPOT(hTotalPot->Integral());
}


// Returns a NuMAtrixSpectrum
NuMatrixSpectrum NuMatrix1D::GetNuMatrixSpectrum(void) const
{
  return NuMatrixSpectrum(*fSpectrum, GetPOT());
}


void NuMatrix1D::Divide(const NuMatrix &by)
{
  // Convert the NuMatrix to a NuMatrix1D
  const NuMatrix1D *other = dynamic_cast<const NuMatrix1D*>(&by);
  if (other == 0) {
    MSG("NuMatrix1D",Msg::kWarning) << "Attempting to divide a NuMatrix1D to something other than a NuMatrix1D" << endl;
  }
  
  this->Divide(*other);
}


void NuMatrix1D::Divide(const NuMatrix1D &by)
{
  fSpectrum->Divide(by.Spectrum());
}


void NuMatrix1D::Divide(const TH1& correction, const Option_t *option)
{
  fSpectrum->Divide(fSpectrum, &correction, 1, 1, option);
}


void NuMatrix1D::Divide(const TGraph &correction)
{
  // copied form NuMatrixSpectrum
  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,""));
  }
  
}


void NuMatrix1D::Multiply(const NuMatrix &by)
{
  // Convert the NuMatrix to a NuMatrix1D
  const NuMatrix1D *other = dynamic_cast<const NuMatrix1D*>(&by);
  if (other == 0) {
    MSG("NuMatrix1D",Msg::kWarning) << "Attempting to multiply a NuMatrix1D to something other than a NuMatrix1D" << endl;
  }
  
  this->Multiply(*other);
}


void NuMatrix1D::Multiply(const NuMatrix1D &by)
{
  fSpectrum->Multiply(by.Spectrum());
}


void NuMatrix1D::Multiply(const TH1& correction, const Option_t *option)
{
  fSpectrum->Multiply(fSpectrum, &correction, 1, 1, option);
}


void NuMatrix1D::Multiply(const TGraph &correction)
{
  // copied form NuMatrixSpectrum
  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,""));
  }
  
}


void NuMatrix1D::Oscillate(const Double_t dm2, const Double_t sn2)
{
  // Copy the simple oscillation function from NuMatrixSpectrum for now
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    Double_t energy = fSpectrum->GetBinCenter(i);
    Double_t oscProb =
          1. - sn2*TMath::
          Power(TMath::Sin(1.27*dm2*(NuUtilities::FDDistance()/Munits::km)
          /energy),2);
    fSpectrum->SetBinContent(i,fSpectrum->GetBinContent(i)*oscProb);
    fSpectrum->SetBinError(i,fSpectrum->GetBinError(i)*oscProb);
  }
}


Double_t NuMatrix1D::StatsLikelihood(const NuMatrix1D* with) const
{
  Double_t like = 0;

  //Aim to minimise -2lnL. That's what I'm returning.
  for (Int_t i=1; i<= fSpectrum->GetNbinsX(); ++i)
  {
    Double_t mnu = fSpectrum->GetBinContent(i);
    Double_t dnu = with->Spectrum()->GetBinContent(i);

    like += NuUtilities::LogLikelihood(dnu,mnu);
  }
  return like;
}



void NuMatrix1D::SetValue(Double_t val)
{
  // No spectrum - don't set
  if (!fSpectrum) {
    MSG("NuMatrix1D",Msg::kWarning) << "Cannot set constant value: No histogram" << endl;
    return;
  }
  
  for (Int_t i = 1; i <= fSpectrum->GetNbinsX(); ++i)
  {
    fSpectrum->SetBinContent(i, val);
  }
}


Int_t NuMatrix1D::Write(const TString &name)
{
  if (fSpectrum) {
    if (name == "")
    {
      return fSpectrum->Write();
    } else {
      return fSpectrum->Write(name);
    }
  } else {
    return 0;
  }
}
  

void NuMatrix1D::Scale(Double_t value)
{
  if (fSpectrum) fSpectrum->Scale(value);
}


void NuMatrix1D::RecoToTrue(const TH2D &correction)
{
  // MSG("NuMatrix1D",Msg::kWarning) << "RecoToTrue not implemented uet for NuMatrix1D" << endl;
  this->YToX(correction);
}


void NuMatrix1D::TrueToReco(const TH2D &correction)
{
  // MSG("NuMatrix1D",Msg::kWarning) << "TrueToReco not implemented uet for NuMatrix1D" << endl;
  this->XToY(correction);
}


void NuMatrix1D::ExtrapolateNDToFD(const TH2D &correction)
{
  // MSG("NuMatrix1D",Msg::kWarning) << "ExtrapolateNDToFD not implemented uet for NuMatrix1D" << endl;
  this->XToY(correction);
}


void NuMatrix1D::XToY(const TH2D &correction)
{
  // This routine taken from NuMatrixspectrum - looks suspect (especially
  // with the differences between the two functions). Do quick validation.
  assert(fSpectrum->GetNbinsX() == correction.GetNbinsX());
  assert(correction.GetNbinsX() == correction.GetNbinsY());
  
  //working space:
  Double_t *val = new Double_t[fSpectrum->GetNbinsX()+1];  
  for(int i=1;i<=fSpectrum->GetNbinsX();i++){
    val[i-1] = 0;
  }

  for(int i=1;i<=fSpectrum->GetNbinsX();i++){ //loop over true
    for(int j=1;j<=fSpectrum->GetNbinsX()+1;j++){ //loop over reco
      val[j-1] += ( fSpectrum->GetBinContent(i) * correction.GetBinContent(i,j));
    }
  }
  for(int i=1;i<=fSpectrum->GetNbinsX()+1;i++) {
    fSpectrum->SetBinContent(i,val[i-1]);
  }
  if (val) {delete[] val; val = 0;}
}

//____________________________________________________________________72
void NuMatrix1D::YToX(const TH2D &correction)
{
  // This routine taken from NuMatrixspectrum - looks suspect (especially
  // with the differences between the two functions). Do quick validation.
  assert(fSpectrum->GetNbinsX() == correction.GetNbinsX());
  assert(correction.GetNbinsX() == correction.GetNbinsY());
  
  //working space:
  Double_t *val = new Double_t[fSpectrum->GetNbinsX()+1];
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) { val[i-1] = 0; }
  
  for(int i=1;i<=fSpectrum->GetNbinsX();i++){
    for(int j=1;j<=correction.GetNbinsY();j++){
      val[j-1] += ( fSpectrum->GetBinContent(i) * correction.GetBinContent(j,i) );
    }
  }
  for(int i=1;i<=fSpectrum->GetNbinsX();i++) {
    fSpectrum->SetBinContent(i,val[i-1]);
  }
  if (val) {delete[] val; val = 0;}
  return;
}

---