NuFCGridPoint.cxx

Go to the documentation of this file.

#include <string>
#include <sstream>
#include <vector>
#include <map>  // For std::pair
#include <stdexcept>


using namespace std;

#include "TH1.h"
#include "TRandom3.h"

#include "NtupleUtils/NuMMRunFC.h"
#include "NtupleUtils/NuMatrixSpectrum.h"

#include "NtupleUtils/NuFCGridPoint.h"
#include "NtupleUtils/NuFCExperiment.h"
#include "NtupleUtils/NuFCExperimentFactory.h"
#include "NtupleUtils/NuFCFitter.h"
#include "NtupleUtils/NuFCConfig.h"

#include "MessageService/MsgService.h"

CVSID("$Id: NuFCGridPoint.cxx,v 1.10 2009/04/15 13:53:41 nickd Exp $");


NuFCGridPoint::NuFCGridPoint(Double_t sin, Double_t dm2, const NuMatrixSpectrum &ndNQData,
  const NuMatrixSpectrum &ndPQData, const NuMMHelperCPT &helper)
  : fsin2(sin),
    fdm2(dm2),
    fHelper(helper),
    ffdPredictionNQ(0),
    ffdPredictionPQ(0),
    fndNQ(ndNQData),
    fndPQ(ndPQData),
    fRandy(0),
    fArchive(false),
    fArchiveFilename("archive.root"),
    fExperimentFactory(0)
{ 
  // Initialise the NuMMParameters object
  // Use default parameters for neutrinos, even though this is bad
  fmmPars.Dm2(2.43e-3);
  fmmPars.Sn2(1.0);

  // Set the neutrino parameters
  fmmPars.Dm2Bar(fdm2);
  fmmPars.Sn2Bar(fsin2);
  // Fix everything except antineutrino fitting
  fmmPars.FixNorm();
  fmmPars.FixShwScale();
  fmmPars.FixNCBackground();
  fmmPars.FixDm2();
  fmmPars.FixSn2();

  // Set up the random number generator
  fRandy = new TRandom3(0);
  
  // Load the configuration
  NuFCConfig &cfg = NuFCConfig::GetConfig();
  // Set up the experiment factory... if we want to

  // We always use the experiment factory now
  string near_events = cfg.GetEventsND();
  string far_events = cfg.GetEventsFD();
  string tau_events = cfg.GetEventsTau();
  fExperimentFactory = new NuFCExperimentFactory(fRandy, near_events, far_events, tau_events);
  
  // Turn off systematics if we've said we don't want them
  if (cfg.Systematics() == false) fExperimentFactory->NoSystematics();
  
  // Are we doing a one-parameter fit?
  if (cfg.FixedSin()) {
    fmmPars.FixSn2Bar();
  }
}


NuFCGridPoint::~NuFCGridPoint()
{
  // Delete the far detector predictions
  if (ffdPredictionNQ) delete ffdPredictionNQ;
  if (ffdPredictionPQ) delete ffdPredictionPQ;
  ffdPredictionNQ = ffdPredictionPQ = 0;
  
  // // Delete the near detector spectra
  // if (fndNQ) delete fndNQ;
  // if (fndPQ) delete fndPQ;
  // fndNQ = fndPQ = 0;

}


// Uses user-generated spectra for the far detector data
void NuFCGridPoint::UsePrediction(const NuMatrixSpectrum &fdNQData,
                                  const NuMatrixSpectrum &fdPQData)
{
  // Delete the far spectra if they exist
  if (ffdPredictionNQ) delete ffdPredictionNQ;
  if (ffdPredictionPQ) delete ffdPredictionPQ;

  // Make a copy of the passed in data spectra
  // this->ffdPredictionNQ = static_cast<TH1*>( fdNQData.Spectrum()->Clone() );
  // this->ffdPredictionPQ = static_cast<TH1*>( fdPQData.Spectrum()->Clone() );
  // ffdPredictionNQ->SetDirectory(0);
  // ffdPredictionPQ->SetDirectory(0);
  ffdPredictionNQ = new NuMatrixSpectrum(fdNQData);
  ffdPredictionPQ = new NuMatrixSpectrum(fdPQData);
}


void NuFCGridPoint::PredictSpectrum(Double_t POT)
{
  // We want to make an extrapolation of the near detector data
  MSG("NuFCGridPoint", Msg::kInfo) << "Generating Extrapolated Near Detector\n"
    << "Using Oscillation Parameters sn2bar: " << fmmPars.Sn2Bar()
    << ", dm2bar: " << fmmPars.Dm2Bar() << endl;
  
  // Give blank prediction spectra, as we won't be checking the far data
  NuMatrixSpectrum blankpot(POT);
  // Make the extrapolator
  NuMMRunFC mmFakeFar(&fHelper, &fndNQ, &fndPQ, &blankpot, &blankpot);
  mmFakeFar.QuietModeOn();
  
  // Now, grab the extrapolations, and pretend they are far detector data
  pair<NuMatrixSpectrum,NuMatrixSpectrum> fakefar = mmFakeFar.MakeFDPred(fmmPars);

  // Tell the user what POT we got out
  MSG("NuFCGridPoint", Msg::kInfo) << "Far extrapolation is at POT "
    << fakefar.second.PoT() << " (" << fakefar.second.Spectrum()->Integral()
    << " events)" << endl;

  // Grab a copy of these Spectra
  // ffdPredictionNQ = static_cast<TH1*>( fakefar.first.Spectrum()->Clone() );
  // ffdPredictionPQ = static_cast<TH1*>( fakefar.second.Spectrum()->Clone() );
  ffdPredictionNQ = new NuMatrixSpectrum(fakefar.first);
  ffdPredictionPQ = new NuMatrixSpectrum(fakefar.second);
}


void NuFCGridPoint::SetRandom(TRandom3 *ranGen)
{
  // Reassign our random number generator, deleting if needed
  if (fRandy) delete fRandy;
  fRandy = ranGen;
}


void NuFCGridPoint::SetSeed(UInt_t seed)
{
  fRandy->SetSeed(seed);
}


const NuMatrixSpectrum &NuFCGridPoint::GetBarSpectrum()
{
  if (!ffdPredictionPQ) {
    // We are backed into a corner... throw an exception as this is
    // a better alternative than a segfault!
    ostringstream err;
    err << "NuFCGridPoint::GetBarSpectrum(): Asking for reference to prediction\n"
        << "    before it is calculated. Try calling NuFCGridPoint::PredictSpectrum\n"
        << "    or NuFCGridPoint::UsePrediction first.";
    throw runtime_error(err.str().c_str());
  }
   return *ffdPredictionPQ;
}

const NuMatrixSpectrum &NuFCGridPoint::GetNuSpectrum()
{
  if (!ffdPredictionPQ) {
    // Also backed into a corner here. Nothing left but an exception!
    ostringstream err;
    err << "NuFCGridPoint::GetNuSpectrum(): Asking for reference to prediction\n"
        << "    before it is calculated. Try calling NuFCGridPoint::PredictSpectrum\n"
        << "    or NuFCGridPoint::UsePrediction first.";
    throw runtime_error(err.str().c_str());
  }
  
  return *ffdPredictionPQ;
}


void NuFCGridPoint::Run(UInt_t number)
{
  for (UInt_t i = 0; i < number; i++)
  {
    // Say we are starting!
    MSG("NuFCGridPoint", Msg::kInfo) << "Running experiment " << i << endl;
    
    // Create a new experiment
    // NuFCExperiment fcexp(i, *ffdPredictionPQ, fRandy);

    // Make a blank extrapolator for predictions
    NuMatrixSpectrum blankpot(ffdPredictionPQ->PoT());
    NuMMRunFC extrapolator(&fHelper, &fndNQ, &fndPQ, &blankpot, &blankpot);
    extrapolator.QuietModeOn();

    // Generte the fake experiment
    fExperimentFactory->GenerateNewExperiment(&extrapolator, fmmPars);
    NuFCExperiment fcexp(i, fExperimentFactory->GetNQSpectrum(), fExperimentFactory->GetPQSpectrum());
    
    // Grab the (possibly shifted) near detector histograms
    NuMatrixSpectrum ndNQ = fExperimentFactory->GetNQSpectrumND();
    NuMatrixSpectrum ndPQ = fExperimentFactory->GetPQSpectrumND();
    
    // Create the fitter
    NuFCFitter fitter(&fcexp, &fHelper, &ndNQ, &ndPQ, fmmPars);
    
    // Do we want to archive every fit?
    if (fArchive) fitter.AlwaysArchive(fArchiveFilename);
    else fitter.ArchiveTo(fArchiveFilename);
    
    // Run the fit
    Double_t fitDeltachi = fitter.Fit(fmmPars);
    
    // Store the delta chi2
    fdeltaChi2.push_back(fitDeltachi);
  }
}

void NuFCGridPoint::SetPars(const NuMMParameters &pars)
{
    fmmPars = pars;    
}

---