NuEZFitter.cxx

Go to the documentation of this file.

#include <string>
#include <map>
#include <iostream>

#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuMMHelperCPT.h"
#include "NtupleUtils/NuMMRunFC.h"
#include "NtupleUtils/NuSystFitter.h"
#include "NtupleUtils/NuMMParameters.h"

#include "NtupleUtils/NuEZFitter.h"

#include "MessageService/MsgService.h"

CVSID("$Id: NuEZFitter.cxx,v 1.4 2009/10/15 16:41:38 ahimmel Exp $");

using namespace std;

void NuEZFitter::Init(void)
{
  fNDDataPQ = fNDDataNQ = fFDDataPQ = fFDDataNQ = 0;
  fHelper = 0; fRun = 0; fFitter = 0; fFitResult = 0;
  
  // Make a dummy fit result
  fFitResult = new NuMMParameters();
  fFitResult->MinuitPass(false);
}

void NuEZFitter::Init_Helper_Near(string Helper, string NeardetData)
{
  // Reading near detector data
  MSG("NuEZFitter",Msg::kInfo) << "Reading Near detector data " << NeardetData << endl;
  fNDDataPQ = new NuMatrixSpectrum(NeardetData, "RecoEnergyPQ_ND");
  fNDDataNQ = new NuMatrixSpectrum(NeardetData, "RecoEnergy_ND");
  
  // Grab the helper file
  MSG("NuEZFitter",Msg::kInfo) << "Extrapolating with helper " << Helper << endl;
  string xsec = "$SRT_PRIVATE_CONTEXT/NtupleUtils/data/xsec_minos_modbyrs4_v3_5_0_mk.root";
  fHelper = new NuMMHelperCPT(Helper.c_str(),xsec);
}
NuEZFitter::NuEZFitter(string Helper, string NeardetData) // :
 //  fNDDataPQ(0), fNDDataNQ(0), fFDDataPQ(0), fFDDataNQ(0),
 //  fHelper(0), fRun(0), fFitter(0), fFitResult(0)
{
  Init();
  Init_Helper_Near(Helper, NeardetData);
  Rebuild();
}

NuEZFitter::NuEZFitter(string Helper, string NeardetData, string FardetData)
{
  Init();
  Init_Helper_Near(Helper, NeardetData);
  UseFDData(FardetData);
}

NuEZFitter::NuEZFitter(const NuEZFitter &o)
{ 
  Init();
  
  if (o.fNDDataPQ) fNDDataPQ = new NuMatrixSpectrum(*o.fNDDataPQ);
  if (o.fNDDataNQ) fNDDataNQ = new NuMatrixSpectrum(*o.fNDDataNQ);
  if (o.fFDDataPQ) fFDDataPQ = new NuMatrixSpectrum(*o.fFDDataPQ);
  if (o.fFDDataNQ) fFDDataNQ = new NuMatrixSpectrum(*o.fFDDataNQ);
  
  if (o.fHelper) fHelper = new NuMMHelperCPT(*o.fHelper);
  if (o.fRun) fRun = new NuMMRunFC(*o.fRun);
  if (o.fFitter) fFitter = new NuSystFitter(*o.fFitter);
  if (o.fFitResult) fFitResult = new NuMMParameters(*o.fFitResult);
}

NuEZFitter::~NuEZFitter()
{
  // Free all of our memory
  if (fNDDataPQ) delete fNDDataPQ; fNDDataPQ = 0;
  if (fNDDataNQ) delete fNDDataNQ; fNDDataNQ = 0;
  
  if (fFDDataPQ) delete fFDDataPQ;  fFDDataPQ = 0;
  if (fFDDataNQ) delete fFDDataNQ;  fFDDataNQ = 0;

  if (fHelper) delete fHelper; fHelper = 0;
  if (fRun) delete fRun; fRun = 0;
  if (fFitter) delete fFitter; fFitter = 0;

  if (fFitResult) delete fFitResult; fFitResult = 0;
// if () delete ;  = 0;

}

//void NuEZFitter::Make_Blank_Fardet(void)
//{
//
//}

void NuEZFitter::Rebuild(void)
{
  // Check the needed plots exist
  if (!(fHelper || fNDDataNQ || fNDDataPQ || fFDDataNQ || fFDDataPQ)) {
    // We don't have all the required plots - just return without doing anything
    MSG("NuEZFitter",Msg::kWarning) << "Rebuild() was called without all the plots in place" << endl;
    return;
  }
  
  // Delete the old objects if necessary
  if (fRun) delete fRun; fRun = 0;
  if (fFitter) delete fFitter; fFitter = 0;
  
  // Build the run object
  // Use a different initialisation if we have no far detector data
  if (!(fFDDataPQ && fFDDataNQ)) {
    fRun = new NuMMRunFC(fHelper, fNDDataNQ, fNDDataPQ, 0);
  } else {
    fRun = new NuMMRunFC(fHelper, fNDDataNQ, fNDDataPQ, fFDDataNQ, fFDDataPQ);
  }

  // Limit the maximum fitting energy
  fRun->SetMaximumEnergy(49.99999);
  
  // And the fitter
  fFitter = new NuSystFitter();
  fFitter->push_back(fRun);
}

void NuEZFitter::MakeFDData(Double_t POT, Double_t dm2bar, Double_t sn2bar, Double_t dm2, Double_t sn2)
{
    // Make the extrapolator
    NuMatrixSpectrum blankpot(POT);
    NuMMRunFC mmFakeFar(fHelper, fNDDataNQ, fNDDataPQ, &blankpot, &blankpot);
    mmFakeFar.QuietModeOn();
    mmFakeFar.SetMaximumEnergy(49.999);
 
    NuMMParameters mmPars;
    mmPars.Dm2(dm2);
    mmPars.Sn2(sn2);
    mmPars.Dm2Bar(dm2bar);
    mmPars.Sn2Bar(sn2bar);
    
    // Clear any existing far detector spectra
    if (fFDDataPQ) delete fFDDataPQ; fFDDataPQ = 0;
    if (fFDDataNQ) delete fFDDataNQ; fFDDataNQ = 0;

    // Grab the extrapolations
    // NuMatrixSpectrum *nqms = 0;
    // NuMatrixSpectrum *pqms = 0;
    mmFakeFar.MakeFDPredNoPair(mmPars, &fFDDataNQ, &fFDDataPQ);
    
    // Rebuild the internal objects
    Rebuild();
}

void NuEZFitter::UseFDData(string FardetData)
{
  // Clear any existing far detector spectra
  if (fFDDataPQ) delete fFDDataPQ; fFDDataPQ = 0;
  if (fFDDataNQ) delete fFDDataNQ; fFDDataNQ = 0;

  MSG("NuEZFitter",Msg::kInfo) << "Reading Far detector data " << FardetData << endl;
  fFDDataPQ = new NuMatrixSpectrum(FardetData, "RecoEnergyPQ_FD");
  fFDDataNQ = new NuMatrixSpectrum(FardetData, "RecoEnergy_FD");
  
  // Rebuild the internal objects
  Rebuild();
}

void NuEZFitter::UseFDData(const NuMatrixSpectrum &fdNQdata, const NuMatrixSpectrum &fdPQdata)
{
  // Clear any existing far detector spectra
  if (fFDDataPQ) delete fFDDataPQ; fFDDataPQ = 0;
  if (fFDDataNQ) delete fFDDataNQ; fFDDataNQ = 0;
  
  fFDDataNQ = new NuMatrixSpectrum(fdNQdata);
  fFDDataPQ = new NuMatrixSpectrum(fdPQdata);
  
  // Now rebuild the fitters
  Rebuild();
}
  
void NuEZFitter::UseFDFakeData(Double_t POT, std::string FardetData, std::string FarDetTauData)
{
  // Clear any existing far detector spectra
  if (fFDDataPQ) delete fFDDataPQ; fFDDataPQ = 0;
  if (fFDDataNQ) delete fFDDataNQ; fFDDataNQ = 0;
  
  // Combine the fake far data and taudata then use as the spectrum
  MSG("NuEZFitter",Msg::kInfo) << "Reading Far detector fake-data " << FardetData << endl;
  MSG("NuEZFitter",Msg::kInfo) << "Reading Far detector fake-Tau data " << FarDetTauData << endl;
  
  NuMatrixSpectrum *fdDataPQ = new NuMatrixSpectrum(FardetData, "RecoEnergyPQ_ND");
  NuMatrixSpectrum *fdDataNQ = new NuMatrixSpectrum(FardetData, "RecoEnergy_ND");
  NuMatrixSpectrum *fdTauDataPQ = new NuMatrixSpectrum(FarDetTauData, "RecoEnergyPQ_ND");
  NuMatrixSpectrum *fdTauDataNQ = new NuMatrixSpectrum(FarDetTauData, "RecoEnergy_ND");

  if (fdDataPQ->PoT() < POT || fdTauDataPQ->PoT() < POT) {
    MSG("NuEZFitter",Msg::kWarning) << "Fake Data files have lower POT than requested - scaling data up" << endl;
  }
  
  fdDataPQ->ScaleToPot(POT);
  fdDataNQ->ScaleToPot(POT);
  fdTauDataPQ->ScaleToPot(POT);
  fdTauDataNQ->ScaleToPot(POT);
  
  fdDataPQ->Add(*fdTauDataPQ);
  fdDataNQ->Add(*fdTauDataNQ);
  
  fFDDataNQ = fdDataNQ;
  fFDDataPQ = fdDataPQ;
  
  // Rebuild the internal objects
  Rebuild();
}

void NuEZFitter::SetPoT(Double_t POT)
{
  // If we don't have far detector spectra,
  if (!fFDDataPQ || !fFDDataNQ) {
    MSG("NuEZFitter",Msg::kWarning) << "No far detector spectra exists. Creating blank spectra" << endl;
    // In case we somehow have one of the two
    if (fFDDataPQ) delete fFDDataPQ; fFDDataPQ = 0;
    if (fFDDataNQ) delete fFDDataPQ; fFDDataPQ = 0;
    
    fFDDataPQ = new NuMatrixSpectrum(POT);
    fFDDataNQ = new NuMatrixSpectrum(POT);
  } else {
    fFDDataPQ->ResetPoT(POT);
    fFDDataNQ->ResetPoT(POT);
  }

  // Far detector spectra have changed... rebuild
  Rebuild();
}

void NuEZFitter::ScaleToPot(Double_t POT)
{
    // If we don't have far detector spectra,
    if (!fFDDataPQ || !fFDDataNQ) {
        MSG("NuEZFitter",Msg::kWarning) << "No far detector spectra exists. Cannot rescale pot." << endl;
        return;
    } else {
        fFDDataPQ->ScaleToPot(POT);
        fFDDataNQ->ScaleToPot(POT);
    }
    
    // Far detector spectra have changed... rebuild
    Rebuild();
}

Double_t NuEZFitter::Chi2(Double_t dm2bar, Double_t sn2bar, Double_t dm2, Double_t sn2) const
{
  NuMMParameters mmPars;
  mmPars.Dm2(dm2);
  mmPars.Sn2(sn2);
  mmPars.Dm2Bar(dm2bar);
  mmPars.Sn2Bar(sn2bar);
  
  // return (*fFitter)(mmPars.VectorParameters());
  return Chi2(mmPars);
}

Double_t NuEZFitter::Chi2(const NuMMParameters &params) const
{
  if (fFitter) {
    return (*fFitter)(params.VectorParameters());
  } else {
    MSG("NuEZFitter",Msg::kError) << "No Fitter exists. Have you added FD data?" << endl;
    return 0;
  }
}

NuMatrixSpectrum NuEZFitter::PredictionNQ(Double_t dm2bar, Double_t sn2bar,
                                        Double_t dm2, Double_t sn2) const
{
  NuMMParameters mmPars;
  mmPars.Dm2(dm2);
  mmPars.Sn2(sn2);
  mmPars.Dm2Bar(dm2bar);
  mmPars.Sn2Bar(sn2bar);
  
  return PredictionNQ(mmPars);
}

NuMatrixSpectrum NuEZFitter::PredictionPQ(Double_t dm2bar, Double_t sn2bar,
                                        Double_t dm2, Double_t sn2) const
{
  NuMMParameters mmPars;
  mmPars.Dm2(dm2);
  mmPars.Sn2(sn2);
  mmPars.Dm2Bar(dm2bar);
  mmPars.Sn2Bar(sn2bar);
  
  return PredictionPQ(mmPars);
}

NuMatrixSpectrum NuEZFitter::PredictionNQ(const NuMMParameters &params) const
{
  // Generate a prediction and return it in a memory-safe way
  NuMatrixSpectrum *pred = fRun->MakeFDPredNuMu(params);
  NuMatrixSpectrum ret(*pred);
  if (pred) delete pred; pred = 0;
  
  return ret;
}

NuMatrixSpectrum NuEZFitter::PredictionPQ(const NuMMParameters &params) const
{
  // Generate a prediction and return it in a memory-safe way
  NuMatrixSpectrum *pred = fRun->MakeFDPredNuBar(params);
  NuMatrixSpectrum ret(*pred);
  if (pred) delete pred; pred = 0;
  
  return ret;
}

Bool_t NuEZFitter::Fit(const NuMMParameters &Initial)
{
  // Can we do a fit?
  if (!fFitter) {
    MSG("NuEZFitter",Msg::kError) << "No Fitter exists. Have you added FD data?" << endl;
    fFitResult = new NuMMParameters(Initial);
    return false;
  }
  
  // Do the fit
  // Bool_t minuitresult = false;
  NuMMParameters mmFit = fFitter->Minimise(Initial);
  // Copy the fit result to our class member
  if (fFitResult) delete fFitResult; fFitResult = 0;
  fFitResult = new NuMMParameters(mmFit);
  
  return fFitResult->MinuitPass();
}
  

---