NuFCFitter.cxx

Go to the documentation of this file.

#include <iostream>
#include <stdexcept>

using namespace std;

#include "TH2D.h"
#include "TFile.h"

#include "NtupleUtils/NuMMParameters.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuFCExperiment.h"
#include "NtupleUtils/NuFCFitter.h"
#include "NtupleUtils/NuMMRunFC.h"
#include "NtupleUtils/NuFCConfig.h"


#include "MessageService/MsgService.h"

CVSID("$Id: NuFCFitter.cxx,v 1.11 2009/06/03 12:51:13 nickd Exp $");


NuFCFitter::NuFCFitter(
    NuFCExperiment *experiment,
    NuMMHelperCPT *helper,
    NuMatrixSpectrum *ndNQData,
    NuMatrixSpectrum *ndPQData,
    const NuMMParameters &trueparams ) 
  : fRun(0),
    ffdNQData(0),
    ffdPQData(0),
    fsystFitter(),
    fBestFit(0),
    fFineGridSearch(false),
    fNegativeDeltaChi(false),
    fHighFitPoint(false),
    fZeroFitPoint(false),
    fArchive(false),
    fArchiveFilename("archive.root"),
    fTrueParams(trueparams)
{
    // Firstly, get copies of the far detector spectrum
  //NuMatrixSpectrum *fdNQ = &experiment->GetSpectrum();
  ffdPQData = new NuMatrixSpectrum(experiment->GetPQSpectrum());
  ffdNQData = new NuMatrixSpectrum(experiment->GetNQSpectrum());
  
  // Create the run!
  fRun = new NuMMRunFC(helper, ndNQData, ndPQData, ffdNQData, ffdPQData);
  fRun->PredictNus(false);
  fRun->SetMaximumEnergy(49.999);
  // fRun->QuietModeOn();
  
  // Now add this run to the fitter
  fsystFitter.push_back(fRun);
  fsystFitter.QuietModeOn();
}


NuFCFitter::~NuFCFitter()
{
  // Delete our run object
  if (fRun) delete fRun;
  fRun = 0;
  
  // Delete our copies of the far detector data
  if (ffdNQData) delete ffdNQData;
  if (ffdPQData) delete ffdPQData;
  ffdNQData = ffdPQData = 0;
  
  // Delete the copy of the best fit point
  if (fBestFit) delete fBestFit;
  fBestFit = 0;
}


Double_t NuFCFitter::Fit(const NuMMParameters &fitParameters)
{
  // Start with a coarse grid search, to seed minuit
  const NuMMParameters mmGrid = CoarseGridSearch(fitParameters);    

  // Now do a minuit fit
  NuMMParameters mmFit = MinuitFit(mmGrid);

  Double_t dchi = DeltaChi(mmFit);

  // If this gave a negative value, rescue it (allow a tiny gap)
  if (dchi < 0) {
    if (dchi < -1e-5) {
      mmFit = RecoverNegativeDeltaChi(mmFit);
      dchi = DeltaChi(mmFit);
    } else {
      // It is below zero, but within error. Set it to zero!
      dchi = 0;
    }
  }

  // Make a copy of the best fit point
  if (fBestFit) delete fBestFit;
  fBestFit = new NuMMParameters(mmFit);

  // Classify the best fit point i.e. is it zero, high etc
  // this is mainly for informational purposes (i.e. the run script)
  ClassifyFitPoint(mmFit);

  // Do we want to archive?
  if (fFineGridSearch || fArchive) { // || fNegativeDeltaChi
    Archive(fArchiveFilename);
  }

  // Do the output
  DoOutput(mmFit, dchi);
  
  return dchi;
}


void NuFCFitter::DoOutput(const NuMMParameters &mmFit, Double_t deltachi )
{
  // Output the fit results
  Double_t fitlike = fsystFitter(mmFit.VectorParameters());
  MSG("NuFCFitter", Msg::kInfo)
    << "    Minimum is dm2bar: " << mmFit.Dm2Bar()
    << " sn2bar: " << mmFit.Sn2Bar() << '\n'
    << "    Fit Likelyhood: " << fitlike << '\n'
    << "    Delta Chi2:     " << deltachi << '\n';
  
  // Did we want to ask the control script to log this?
  if (fFineGridSearch || fNegativeDeltaChi)
    MSG("NuFCFitter", Msg::kInfo) << "    Please Log this experiment." << endl;
  
}


NuMMParameters NuFCFitter::MinuitFit(NuMMParameters mmStart)
{
  // If it is trying to start us at zero, reset to the search minimum
  // fTrueParams.LowerDm2BarLimit()) {
  if (mmStart.Sn2Bar() == 0) {
    MSG("NuFCFitter", Msg::kDebug) << "    Starting minuit at no-oscillation; shifting to lower limits" << endl;
    mmStart.Sn2Bar(fTrueParams.LowerDm2BarLimit());
    mmStart.Dm2Bar(fTrueParams.LowerSn2BarLimit());
  }
  
  // Fix the sin, if specified
  NuFCConfig &cfg = NuFCConfig::GetConfig();
  if (cfg.FixedSin()) {
    mmStart.FixSn2Bar();
    mmStart.Sn2Bar(fTrueParams.Sn2Bar());
  }
  
  // Now run minuit itself
  NuMMParameters mmFit = fsystFitter.Minimise(mmStart);
  
  // Did minuit succeed?
  if (mmFit.MinuitPass())
  {
    // Tell them what minuit found
    MSG("NuFCFitter", Msg::kInfo)
      << "    Minuit found: dm2: " << mmFit.Dm2()
      << " sn2: " << mmFit.Sn2()
      << " dm2bar: " << mmFit.Dm2Bar()
      << " sn2bar: " << mmFit.Sn2Bar()
      << '\n';
  } else {
      // Minuit failed to fit. Use a more intensive grid search,
      // starting on the minuit starting point (that we get from
      // a grid search)
      MSG("NuFCFitter", Msg::kInfo) << "    Minuit fit failed." << endl;
      MSG("NuFCFitter", Msg::kInfo) << "    Searching fine grid around sn2: " << mmStart.Sn2Bar()
                                    << ", dm2: " << mmStart.Dm2Bar() << endl;
      NuMMParameters mmFineGrid = FineGridSearch(mmStart);
      fFineGridSearch = true;
      
      // Test the fine grid search against the result from minuit = it could be possible
      // that minuit found a 'better' point, even if it failed
      Double_t minuit_chi = fsystFitter(mmFit.VectorParameters());
      Double_t finegrid_chi = fsystFitter(mmFineGrid.VectorParameters());
      
      if (minuit_chi > finegrid_chi) {
        // The new, fine grid point is better
        mmFit = mmFineGrid;
      } else {
        // The minuit fail point is actually better. Use that.
        MSG("NuFCFitter", Msg::kInfo) << "    Failed Minuit fit is better than fine grid, using." << endl;
      }
      
  }
  
  return mmFit;
}


NuMMParameters NuFCFitter::CoarseGridSearch(NuMMParameters mmGrid)
{
  // The width and height of the grid
  Int_t gridWidth = 5;
  Double_t sinMin = 0.2;
  Double_t sinMax = 1.0;
  
  const Int_t gridHeight = 15;
  
  // Check if we are doing one-parameter fit, and if so constrain the grid search
  // to one-dimension
  NuFCConfig &cfg = NuFCConfig::GetConfig();
  if (cfg.FixedSin()) {
    gridWidth = 1;
    sinMin = sinMax = fTrueParams.Sn2Bar();
  }

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, sinMin, sinMax, gridHeight, 2e-3, 30e-3);
  mmGrid = GridSearch(likesurf, mmGrid);
  
  // Now do a log search
  vector<Double_t> bins;
  for (Int_t i = -15; i <= -1; i++) {
    bins.push_back(pow(10, (Double_t)i/10.0));
  }
  Int_t bincount = bins.size()-1;
  TH2D logsurf("logsurf", "logsurf", gridWidth, sinMin, sinMax, bincount, &(bins.front()));
  mmGrid = GridSearch(logsurf, mmGrid, &mmGrid);

  // Output the minimum
  MSG("NuFCFitter", Msg::kInfo) << "    Coarse grid search found sn2: "
    << mmGrid.Sn2Bar() << ", dm: " << mmGrid.Dm2Bar() << '\n';
  
  return mmGrid;
}


NuMMParameters NuFCFitter::FineGridSearch(NuMMParameters mmStart)
{
  // Calculate the location of the fine grid search
  Double_t gridL, gridR, gridT, gridB;
  gridL = gridR = gridT = gridB = -1.0;
  
  // The width and height of the grid
  Int_t gridWidth = 19;
  const Int_t gridHeight = 19;
  
  // Configuration
  NuFCConfig &cfg = NuFCConfig::GetConfig();
  
  // change behavior if we want a one-parameter fit
  if (!cfg.FixedSin()) {
      // Go +- 0.2 in sin space
    gridL = mmStart.Sn2Bar() - 0.2;
    gridR = gridL + 0.4;
    // Now limit it to the physical realm
    if (gridL <= fTrueParams.LowerSn2BarLimit()) gridL = fTrueParams.LowerSn2BarLimit();
    if (gridR > 1.0) gridR = 1.0;
  } else {
    // One-parameter fit
    gridL = gridR = fTrueParams.Sn2Bar();
    gridWidth = 1;
  }

  // And go +-2e-3 in mass space
  gridB = mmStart.Dm2Bar() - 2e-3;
  gridT = gridB + 4e-3;
  // We only really need to worry about lower limit here
  if (gridB <= fTrueParams.LowerDm2BarLimit()) gridB = fTrueParams.LowerSn2BarLimit();

  MSG("NuFCFitter", Msg::kInfo) 
    << "    Doing fine grid search on sin = [ "
    << gridL << ", " << gridR << " ]\n"
    << "                              dm  = [ "
    << gridB << ", " << gridT << " ]\n";

  // Ensure that we haven't done anything stupid like set low above high
  if (gridL > gridR || gridT <= gridB) {
//      cout << "Error: Negative range in fine grid search" << endl;
      throw runtime_error("Negative range in fine grid search");
  }
  
  // Make the surface
  
  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
    gridHeight, gridB, gridT);
    
  // And now do a grid search on this
  NuMMParameters mmGrid = GridSearch(likesurf, mmStart);
  
  // Temporarily, archive any fine grid surface
  // TFile tf("archives/finegrid.root", "UPDATE");
  //   likesurf.Write();
  //   tf.Close();
  
  // Output the minimum
  MSG("NuFCFitter", Msg::kInfo) << "    Fine grid search found sn2: "
    << mmGrid.Sn2Bar() << ", dm: " << mmGrid.Dm2Bar() << '\n';
  
  
  return mmGrid;
}


NuMMParameters NuFCFitter::GridSearch(TH2D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref)
{
    Double_t likelihood = fsystFitter.FillLikelihoodSurfaceBar(likesurf, mmGrid);
    
    // Find the point on this surface that is minimum
    Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
    likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
    Double_t minsin = likesurf.GetXaxis()->GetBinCenter(minbinX);
    Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);

    // Access the configuration
    NuFCConfig &cfg = NuFCConfig::GetConfig();

    // Check the zero point, to see if it is lower than anything else on the grid
    mmGrid.Dm2Bar(0);
    // Different 'zero' point if one parameter
    if (cfg.FixedSin()) {
      mmGrid.Sn2Bar(fTrueParams.Sn2Bar());
    } else {
      mmGrid.Sn2Bar(0);
    }
    
    Double_t zeropoint = fsystFitter(mmGrid.VectorParameters());
    if (zeropoint < likelihood) {
        // Then the no-oscillation case is the minimum
        minsin = 0.0;
        minmass = 0.0;
        
        if (cfg.FixedSin()) {
          minsin = fTrueParams.Sn2Bar();
        }
    }
    
    // If we have been given a reference point, check to see if this is lower
    if (Ref) {
      Double_t refpoint = fsystFitter(Ref->VectorParameters());
      if (refpoint < likelihood) {
        MSG("NuFCFitter", Msg::kDebug) << "    Reference point is lower than grid search" << endl;
        minsin = Ref->Sn2Bar();
        minmass = Ref->Dm2Bar();
      }
    }
    
    // Set mmPars to the newly found point, then return it
    mmGrid.Sn2Bar(minsin);
    mmGrid.Dm2Bar(minmass);
    
    return mmGrid;
}
  

Double_t NuFCFitter::DeltaChi(const NuMMParameters& bestfit)
{
    // Calculate the two likelihood points
    Double_t fitlike = fsystFitter(bestfit.VectorParameters());

    Double_t truelikelihood = 0;  

    // cout << "bfdm: " << bestfit.Dm2Bar() << ", limit = " << fTrueParams.LowerDm2BarLimit() << endl;

    // Don't allow the parameters to go below constrained - this can
    // give us an artificially negative delta chi2
    if ( (fTrueParams.Dm2Bar() < fTrueParams.LowerDm2BarLimit()) ||
         (fTrueParams.Sn2Bar() < fTrueParams.LowerSn2BarLimit()) ) {  
        MSG("NuFCFitter", Msg::kDebug) << "    True point is lower than fit is allowed." << endl;
        NuMMParameters pc(bestfit);
        pc.Dm2Bar(fTrueParams.LowerDm2BarLimit());
        pc.Sn2Bar(fTrueParams.LowerSn2BarLimit());
        truelikelihood = fsystFitter(pc.VectorParameters());
    } else {
        truelikelihood = fsystFitter(fTrueParams.VectorParameters());
    }
    

    // fmmPars.Dm2Bar(fmmPars.LowerDm2BarLimit());
    //   MAXMSG("NuFCGriPoint", 5, Msg::kInfo) << "Moving true grid point "
    //   fmmPars.Sn2Bar(fmmPars.LowerSn2BarLimit());
    //   
    // Double_t truelikelihood = fsystFitter(fTrueParams.VectorParameters());

    // Return the difference
    return truelikelihood - fitlike;
}

// We have gotten a negative delta chi2. Recover this!    

NuMMParameters NuFCFitter::RecoverNegativeDeltaChi (const NuMMParameters& bestfit)
{
  // Set the failure flag
  fNegativeDeltaChi = true;

  // Recalculate the likelihoods
  Double_t truelike = fsystFitter(fTrueParams.VectorParameters());
  Double_t oldlike = fsystFitter(bestfit.VectorParameters());;

  // Make a copy of the true parameters, to ensure that we do a 1D fit if specified
  // (this is assuming we cannot guarantee that the true parameters has the parameter fixed)
  NuMMParameters tpcopy(fTrueParams);
  NuFCConfig &cfg = NuFCConfig::GetConfig();
  if (cfg.FixedSin()) {
    tpcopy.FixSn2Bar();
  }
  // NuMMParameters newFit = fsystFitter.Minimise(fTrueParams);
  NuMMParameters newFit = fsystFitter.Minimise(tpcopy);
  
  MSG("NuFCFitter", Msg::kInfo)
    << "    Negative delta chi2 likelihood (" << truelike - oldlike
    << ") found. Refitting from true (" << truelike << ")" << endl;
  
  
  // If this passed, and is lower than our old likelihood, then use it
  if (newFit.MinuitPass())
  {
      MSG("NuFCFitter", Msg::kInfo) << "    Refit passed." << endl;
      Double_t newlike = fsystFitter(newFit.VectorParameters());
      if (newlike > oldlike) {
          MSG("NuFCFitter", Msg::kInfo)
            << "    New fit point ( dm: " << newFit.Dm2Bar() << ", sn: " << newFit.Sn2Bar() << " )\n"
            << "    at " << newlike << " is higher than before ( " << oldlike << " )! strange!\n"
            << "    Returning true fit point, as is lowest we have seen."
            << endl;
          newFit = fTrueParams;
      }
  } else {
      MSG("NuFCFitter", Msg::kError)
        << "    Unrecoverable error: Refit on negative deltachi failed."
        << endl;
      throw logic_error("NuFCFitter Refit on negative deltachi failed.");
  }

  return newFit;
}


void NuFCFitter::ClassifyFitPoint(const NuMMParameters& mmFit)
{
  // cout << "    Checking against minimum dm: " << fTrueParams.LowerDm2BarLimit() << ", sn: " << fTrueParams.LowerSn2BarLimit() << endl;
 
  // Check to see if the fit is to 'zero'
  if (mmFit.Sn2Bar() <= fTrueParams.LowerSn2BarLimit() + 1e-4 || 
      mmFit.Dm2Bar() <= fTrueParams.LowerDm2BarLimit() + 1e-6 )
  {
      MSG("NuFCFitter", Msg::kInfo)
        << "    Best fit is to no oscillation case\n";
      fZeroFitPoint = true;
  }
    
  // Now check if the best fit is maxing out
  // This message is a little out of date, but will do for now
  if (mmFit.Dm2Bar() >= 30e-3 - 0.25e-4) {
    MSG("NuFCFitter", Msg::kInfo)
     << "    Best fit is maxing out on dm2bar ("
     << mmFit.Dm2Bar() << ")" << endl;
    fHighFitPoint = true;
  }
}
  

void NuFCFitter::Archive(std::string filename)
{
  MSG("NuFCFitter",Msg::kInfo) << "Disabled all archiving." << endl;
  return;

  if (!ffdPQData || !ffdNQData) {
    MSG("NUFCFitter",Msg::kError) << "No far detector spectrum: Cannot archive" << endl;
    return;
  }
  if (!fBestFit) {
    MSG("NuFCFitter", Msg::kError)
      << "Error: Trying to archive before fit" << endl;
    return;
  }
  
  MSG("NuFCFitter", Msg::kInfo)
    << "    Archiving to file " << filename << "..." << endl;

  // The actual experiment spectrum
//  TH1D fdSpec(*(ffdPQData->Spectrum()));
  NuMatrixSpectrum fdSpec(*ffdPQData);
//   // fdSpec.RebinTo4GeV();
//   // Grab the name out of this
  string basename = fdSpec.Spectrum()->GetName();
// 
// // No longer do the likelihood surface - it takes too long. Instead just write out the spectrum,
// // and we can regenerate the surface later if we wish
//   
  // Now we want the likelihood surfaces. First do the linear portion
  string surfname = basename + "_surf";
  TH2D likesurface(surfname.c_str(), surfname.c_str(), 50, 0, 1, 100, 0, 20e-3);
//  TH2D likesurface(surfname.c_str(), surfname.c_str(), 100, 0.9, 1, 100, 0.002, 0.003);
//  Double_t minlin = fsystFitter.FillLikelihoodSurfaceBar(likesurface, *fBestFit);
//   
//   // Now do the logarithmic plot
//   vector<Double_t> bins;
//   for (Int_t i = -204; i <= 0; i++) {
//     bins.push_back(pow(10, (Double_t)i/100.0));
//   }
//   Int_t bincount = bins.size()-1;
//   string logname = surfname + "_log";
//   TH2D logsurf(logname.c_str(), logname.c_str(), 50, 0, 1, bincount, &(bins.front()));
//   Double_t minlog = fsystFitter.FillLikelihoodSurfaceBar(logsurf, *fBestFit);
//   
//   // Now, shift both histograms by the minimum
//   Double_t surfmin = minlin < minlog ? minlin : minlog;
//   ShiftHistogram(likesurface, -surfmin);
//   ShiftHistogram(logsurf, -surfmin);
  
  // Now, write these out!
  TFile tf(filename.c_str(), "UPDATE");
    fdSpec.Spectrum()->Write();
    likesurface.Write();
    // logsurf.Write();
  tf.Close();
}
  

void NuFCFitter::ShiftHistogram(TH2D &surface, Double_t shift)
{
    for (Int_t x = 1; x <= surface.GetNbinsX(); x++)
    {
        for (Int_t y = 1; y <= surface.GetNbinsY(); y++)
        {
            Double_t contents = surface.GetBinContent(x, y);
            surface.SetBinContent(x, y, contents + shift);
        }
    }
}

---