NuMatrixCPT.cxx

Go to the documentation of this file.

#include "NuMatrix.h"
#include "NuMatrixCPT.h"

#include "MessageService/MsgService.h"

#include "TMath.h"

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

ClassImp(NuMatrixCPT)


// Constructors

NuMatrixCPT::NuMatrixCPT() : NuMatrix(), fPQMatrix(0), fNQMatrix(0)
{
  
}

NuMatrixCPT::NuMatrixCPT(Double_t POT) : NuMatrix(POT), fPQMatrix(0), fNQMatrix(0)
{
  
}

NuMatrixCPT::NuMatrixCPT(const NuMatrixCPT& orig) : NuMatrix(orig), fPQMatrix(0), fNQMatrix(0)
{
  // Copy the histograms over
  if (orig.fPQMatrix) fPQMatrix = orig.fPQMatrix->Copy();
  if (orig.fNQMatrix) fNQMatrix = orig.fNQMatrix->Copy();
}

NuMatrixCPT &NuMatrixCPT::operator=(const NuMatrixCPT& orig)
{
  // Copy the base class stuff
  this->NuMatrix::operator=(orig);
  
  // If our spectra exist, delete them
  if (fPQMatrix) delete fPQMatrix; fPQMatrix = 0;
  if (fNQMatrix) delete fNQMatrix; fNQMatrix = 0;
  
  // Now, copy over the matrix data from the other object
  if (orig.fPQMatrix) fPQMatrix = orig.fPQMatrix->Copy();
  if (orig.fNQMatrix) fNQMatrix = orig.fNQMatrix->Copy();
  
  return *this;
}

// Load from two seperate NuMatrix objects
NuMatrixCPT::NuMatrixCPT(const NuMatrix* NQ, const NuMatrix* PQ)
  : NuMatrix(),
    fPQMatrix(0),
    fNQMatrix(0)
{
  // Vector to count POT of valid histograms
  vector<Double_t> pots;
  
  if (NQ) {
    fNQMatrix = NQ->Copy();
    pots.push_back(fNQMatrix->GetPOT());
  }
  if (PQ) {
    fPQMatrix = PQ->Copy();
    pots.push_back(fPQMatrix->GetPOT());
  }
  
  // Check that every POT is the same, and display a warning otherwise
  // This is because I can't imagine a legitimate reason to have this
  // scenario, at the moment.
  // Now, how do we set the POT? Check they are both the same, for now.
  // For 'same' define POT to 1e-5. Surely nobody makes values that different?
  if (fPQMatrix && fNQMatrix) {
    if (TMath::Abs(NQ->GetPOT() - PQ->GetPOT()) > 1e-5) {
      MSG("NuMatrixCPT",Msg::kWarning) << "POT count for pre-initialised NuMatrix objects differ." << endl;
      MSG("NuMatrixCPT",Msg::kWarning) << "  PQ: " << fPQMatrix->GetPOT() << endl;
      MSG("NuMatrixCPT",Msg::kWarning) << "  NQ: " << fNQMatrix->GetPOT() << endl;
      MSG("NuMatrixCPT",Msg::kWarning) << "Storing POT of NQ in object" << endl;
    }
  }
  
  // Set out POT to the first one, recorded, if we have any.
  if (pots.size()) ResetPOT(pots[0]);
  
}


// Load histograms from file!
NuMatrixCPT::NuMatrixCPT( const TString &filename,
                          const TString &NQhistname,
                          const TString &PQhistname,
                          const TString &pothistname )
                        : NuMatrix(),
                          fPQMatrix(0),
                          fNQMatrix(0)
{
  // Create the two histograms
  fNQMatrix = new NuMatrix1D(filename, NQhistname, pothistname);
  fPQMatrix = new NuMatrix1D(filename, PQhistname, pothistname);
  
  // Now, Grab the POT out of each
  Double_t potPQ = 0.0, potNQ = 0.0;
  potPQ = fPQMatrix->GetPOT();
  potNQ = fNQMatrix->GetPOT();
  
  if (potPQ != potNQ)
  {
    MSG("NuMatrixCPT",Msg::kWarning) << "POT read from file "
      << filename << " does not seem to match between histograms."
      << endl;
  } else {
    this->ResetPOT(potNQ);
  }
  
}

NuMatrixCPT::~NuMatrixCPT()
{
  if (fPQMatrix) delete fPQMatrix; fPQMatrix = 0;
  if (fNQMatrix) delete fNQMatrix; fNQMatrix = 0;
}


// Normal Functions ////////////////////////////////////////////////////

Double_t NuMatrixCPT::StatsLikelihood(const NuMatrix* with) const
{
  Double_t like = 0.0;
  if (fPQMatrix) like += fPQMatrix->StatsLikelihood(with);
  if (fNQMatrix) like += fNQMatrix->StatsLikelihood(with);
  
  return like;
}


void NuMatrixCPT::ScaleToPOT(Double_t new_POT)
{
  if (fPQMatrix) fPQMatrix->ScaleToPOT(new_POT);
  if (fNQMatrix) fNQMatrix->ScaleToPOT(new_POT);
  this->ResetPOT(new_POT);
}


void NuMatrixCPT::Scale(Double_t value)
{
  if (fPQMatrix) fPQMatrix->Scale(value);
  if (fNQMatrix) fNQMatrix->Scale(value);
}


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


void NuMatrixCPT::Divide(const NuMatrixCPT &correction)
{
  if (fPQMatrix && correction.PQ()) fPQMatrix->Divide(*correction.PQ());
  if (fNQMatrix && correction.NQ()) fNQMatrix->Divide(*correction.NQ());
}


void NuMatrixCPT::Divide(const TH1& correction, const Option_t *option)
{
  if (fPQMatrix) fPQMatrix->Divide(correction, option);
  if (fNQMatrix) fNQMatrix->Divide(correction, option);
}


void NuMatrixCPT::Divide(const TGraph &correction)
{
  if (fPQMatrix) fPQMatrix->Divide(correction);
  if (fNQMatrix) fNQMatrix->Divide(correction);
}


void NuMatrixCPT::Divide(const std::pair<const TH1D*,const TH1D*> &corr)
{
  if (fNQMatrix && corr.first) fNQMatrix->Divide(*corr.first);
  if (fPQMatrix && corr.second) fPQMatrix->Divide(*corr.second);
}


void NuMatrixCPT::Divide(const std::pair<const TGraph*,const TGraph*> &corr)
{
  if (fNQMatrix && corr.first) fNQMatrix->Divide(*corr.first);
  if (fPQMatrix && corr.second) fPQMatrix->Divide(*corr.second);
}


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


void NuMatrixCPT::Multiply(const NuMatrixCPT &correction)
{
  if (fPQMatrix && correction.PQ()) fPQMatrix->Multiply(*correction.PQ());
  if (fNQMatrix && correction.NQ()) fNQMatrix->Multiply(*correction.NQ());  
}


void NuMatrixCPT::Multiply(const TH1& correction, const Option_t *option)
{
  if (fPQMatrix) fPQMatrix->Multiply(correction, option);
  if (fNQMatrix) fNQMatrix->Multiply(correction, option);
}


void NuMatrixCPT::Multiply(const TGraph &correction)
{
  if (fPQMatrix) fPQMatrix->Multiply(correction);
  if (fNQMatrix) fNQMatrix->Multiply(correction);
}

void NuMatrixCPT::Multiply(const std::pair<const TH1D*,const TH1D*> &corr)
{
  if (fNQMatrix && corr.first) fNQMatrix->Multiply(*corr.first);
  if (fPQMatrix && corr.second) fPQMatrix->Multiply(*corr.second);
}


void NuMatrixCPT::Multiply(const std::pair<const TGraph*,const TGraph*> &corr)
{
  if (fNQMatrix && corr.first) fNQMatrix->Multiply(*corr.first);
  if (fPQMatrix && corr.second) fPQMatrix->Multiply(*corr.second);  
}


Int_t NuMatrixCPT::Write(const TString &name)
{
  MSG("NuMatrixCPT",Msg::kWarning) << "Histogram writing properly implemented: Only uses PQ/NQ suffixes." << endl;
  
  Int_t bytetotal = 0;
  
  // If we haven't been given a name, amend the names of the histograms
  if (name=="")
  {
    // Don't allow the user to write with no name - this is because we
    // cannot (currently) access the histograms name in a Storage-independant way
    MSG("NuMatrixCPT",Msg::kWarning) << "Must specify a name when writing. Current limitation of NuMatrix framework." << endl;

    // TString hname;
    // 
    // if (fPQMatrix) {
    //   hname = fPQMatrix->Spectrum()->GetName();
    //   bytetotal += fPQMatrix->Write(hname + "_PQ");
    // }
    // 
    // if (fNQMatrix) {
    //   hname = fNQMatrix->Spectrum()->GetName();
    //   bytetotal += fNQMatrix->Write(name + "_NQ");
    // }
    
  } else {
    if (fPQMatrix) bytetotal += fPQMatrix->Write(name + "_PQ");
    if (fNQMatrix) bytetotal += fNQMatrix->Write(name + "_NQ");
  }
  
  return bytetotal;
}


void NuMatrixCPT::ResetPOT(Double_t new_POT)
{
  fPot = new_POT;
  if (fPQMatrix) fPQMatrix->ResetPOT(new_POT);
  if (fNQMatrix) fNQMatrix->ResetPOT(new_POT);
}
    

void NuMatrixCPT::Draw(Option_t*)
{
  MSG("NuMatrixCPT",Msg::kInfo) << "Cannot draw CPT matrix directly at this time. Draw NQ() and PQ() independently." << endl;
  
}
  

void NuMatrixCPT::Oscillate(const NuMMParameters &p)
{
  Oscillate(p.Dm2(), p.Sn2(),p.Dm2Bar(),p.Sn2Bar());
  
}

void NuMatrixCPT::Oscillate(const Double_t dm2, const Double_t sn2, const Double_t dm2bar, const Double_t sn2bar)
{
  if (fPQMatrix) fPQMatrix->Oscillate(dm2bar, sn2bar);
  if (fNQMatrix) fNQMatrix->Oscillate(dm2, sn2);
}


void NuMatrixCPT::SetValue(Double_t val)
{
  if (fPQMatrix) fPQMatrix->SetValue(val);
  if (fNQMatrix) fNQMatrix->SetValue(val);
}
  

void NuMatrixCPT::RecoToTrue(const TH2D &)
{
  MSG("NuMatrixCPT",Msg::kError) << "Trying to Reco->True a CPT matrix with only one histogram!" << endl;
}


void NuMatrixCPT::TrueToReco(const TH2D &)
{
  MSG("NuMatrixCPT",Msg::kError) << "Trying to True-Reco a CPT matrix with only one histogram!" << endl;  
}


void NuMatrixCPT::ExtrapolateNDToFD(const TH2D &)
{
    MSG("NuMatrixCPT",Msg::kError) << "Trying to extrapolate a CPT matrix with only one histogram!" << endl;
}


void NuMatrixCPT::RecoToTrue(const std::pair<const TH2D*,const TH2D*> &correction)
{
  if (fNQMatrix) fNQMatrix->RecoToTrue(*correction.first);
  if (fPQMatrix) fPQMatrix->RecoToTrue(*correction.second);
}


void NuMatrixCPT::TrueToReco(const std::pair<const TH2D*,const TH2D*> &correction)
{
  if (fNQMatrix) fNQMatrix->TrueToReco(*correction.first);
  if (fPQMatrix) fPQMatrix->TrueToReco(*correction.second);
}


void NuMatrixCPT::ExtrapolateNDToFD(const std::pair<const TH2D*,const TH2D*> &beammatrix)
{
  if (fNQMatrix) fNQMatrix->ExtrapolateNDToFD(*beammatrix.first);
  if (fPQMatrix) fPQMatrix->ExtrapolateNDToFD(*beammatrix.second);
}

---