NuMatrix2D.cxx

Go to the documentation of this file.

// NuMatrix2D.cxx
// A slightly less abstract, 2D data storage class

#include "NuMatrix2D.h"
#include "NuUtilities.h"

#include "MessageService/MsgService.h"

#include "TH2D.h"

CVSID("$Id: NuMatrix2D.cxx,v 1.3 2009/11/30 13:03:10 nickd Exp $");

ClassImp(NuMatrix2D)

void Inst2d(void)
{
  NuMatrix2D a;
}
// Constructors

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

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

// The spectrum-only constructor
NuMatrix2D::NuMatrix2D(const TH2D& spectrum) : NuMatrix(0), fSpectrum(0)
{
  fSpectrum = new TH2D(spectrum);
}

// The spectrum *and* POT constructor
NuMatrix2D::NuMatrix2D(const TH2D& spectrum, Double_t POT) : NuMatrix(POT), fSpectrum(0)
{
  fSpectrum = new TH2D(spectrum);
}

// The copy-constructor
NuMatrix2D::NuMatrix2D(const NuMatrix2D& original) : NuMatrix(original), fSpectrum(0)
{
  fSpectrum = new TH2D(*original.fSpectrum);
}

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


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

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

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


Double_t NuMatrix2D::StatsLikelihood(const NuMatrix2D* with) const
{
  // Do the same as the 1D, but over all bins
  Double_t like = 0.0;
  
  //Aim to minimise -2lnL. That's what I'm returning.
  for (Int_t i=1; i<= fSpectrum->GetNbinsX(); ++i)
  {
    for (Int_t j=1; j <= fSpectrum->GetNbinsY(); ++j)
    {
      //Bizarre limits because root histograms are silly
      Double_t mnu = fSpectrum->GetBinContent(i, j);
      Double_t dnu = with->Spectrum()->GetBinContent(i, j);

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


void NuMatrix2D::ScaleToPOT(Double_t new_POT)
{
  // If the old POT is zero, or the spectrum does not exist
  if (GetPOT() == 0.0 || fSpectrum == 0) {
    // Simply set the new POT (avoids infinite scaling)
    ResetPOT(new_POT);
  } else {  
    // The spectrum must exist - calculate the scale factor
    Double_t scale_factor = new_POT / GetPOT();
    fSpectrum->Scale(scale_factor);
  }
}


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


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


void NuMatrix2D::Divide(const NuMatrix2D &correction)
{
  fSpectrum->Divide(correction.Spectrum());
}


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


void NuMatrix2D::Divide(const TGraph &)
{
  MSG("NuMatrix2D",Msg::kWarning)
    << "2D histograms cannot be corrected with graphs without further information"
    << endl;
}



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


void NuMatrix2D::Multiply(const NuMatrix2D &correction)
{
  fSpectrum->Multiply(correction.Spectrum());
}


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


void NuMatrix2D::Multiply(const TGraph &)
{
  MSG("NuMatrix2D",Msg::kWarning)
    << "2D histograms cannot be corrected with graphs without further information"
    << endl;
}


// Matrix method functions
void NuMatrix2D::RecoToTrue(const TH2D &)
{
  MSG("NuMatrix2D",Msg::kWarning) << "Base NuMatrix2D cannot Reco->True" << endl;
}


void NuMatrix2D::TrueToReco(const TH2D &)
{
  MSG("NuMatrix2D",Msg::kWarning) << "Base NuMatrix2D cannot True->Reco" << endl;  
}


void NuMatrix2D::ExtrapolateNDToFD(const TH2D &)
{
  MSG("NuMatrix2D",Msg::kWarning) << "Base NuMatrix2D cannot be extrapolated" << endl;
}


void NuMatrix2D::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)
  {
    for (Int_t j = 1; j <= fSpectrum->GetNbinsY(); ++j)
    {
      fSpectrum->SetBinContent(i, j, val);
    }
  }
}


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


void NuMatrix2D::Oscillate(const Double_t , const Double_t )
{
  MSG("NuMatrix2D",Msg::kWarning) << "Base NuMatrix2D cannot be oscillated" << endl;  
}
  

void NuMatrix2D::Draw(Option_t* option)
{
  if(fSpectrum) {
    fSpectrum->Draw(option);
  } else {
    MSG("NuMatrix2D",Msg::kWarning) << "No spectrum exists: Cannot draw" << endl;
  }
}
  

---