NuMMRunFC.cxx

Go to the documentation of this file.

#include "TCanvas.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuHistInterpolator.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "NtupleUtils/NuMMHelperCPT.h"
#include "NtupleUtils/NuMMParameters.h"
#include "NtupleUtils/NuMMRunFC.h"

ClassImp(NuMMRunFC)

CVSID("$Id: NuMMRunFC.cxx,v 1.22 2009/11/02 21:32:37 hartnell Exp $");

using namespace Sample;

//____________________________________________________________________72
NuMMRunFC::NuMMRunFC(NuMMHelperCPT* helper,
                     NuMatrixSpectrum* ndNuData,
                     NuMatrixSpectrum* ndBarData,
                     double pot)
               : NuMMRunNuBar(),
                 fPredictNeutrinos(true),
                 fPredictAntiNeutrinos(true),
                 fDoneWsNeutrinos(false),
                 fDoneWsAntineutrinos(false)
{
  TString name = "hError";
  name += (counter++);
  TH1D hBlank(name, "Something Went Wrong", 10, 0, 30);
  
  fndNuData = ndNuData;
  fndBarData = ndBarData;
  fHelper = helper;
  
  ffdNuData = new NuMatrixSpectrum(hBlank, pot);
  ffdBarData = new NuMatrixSpectrum(hBlank, pot);

  
  // Blank out all the Precache variables
  PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
  PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
  PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
  PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;
  
  // We should now be completely capable of calculating the
  // halfway point of the prediction. Do so now.
  CacheExtrapolation();
  
}

//____________________________________________________________________72
NuMMRunFC::NuMMRunFC(NuMMHelperCPT* helper,
                       NuMatrixSpectrum* ndNuData,
                       NuMatrixSpectrum* ndBarData,
                       NuMatrixSpectrum* fdNuData,
                       NuMatrixSpectrum* fdBarData)
               : NuMMRunNuBar(),
                 fPredictNeutrinos(true),
                 fPredictAntiNeutrinos(true),
                 fDoneWsNeutrinos(false),
                 fDoneWsAntineutrinos(false)
{
  fndNuData = ndNuData;
  fndBarData = ndBarData;
  ffdNuData = fdNuData;
  ffdBarData = fdBarData;
  fHelper = helper;
  
  // Blank out all the Precache variables
  PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
  PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
  PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
  PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;
  
  // We should now be completely capable of calculating the
  // halfway point of the prediction. Do so now.
  CacheExtrapolation();
  
}


//____________________________________________________________________72
// The copy contructor
NuMMRunFC::NuMMRunFC(const NuMMRunFC &old) : NuMMRunNuBar()
{
  fndNuData = old.fndNuData;
  fndBarData = old.fndBarData;
  ffdNuData = old.ffdNuData;
  ffdBarData = old.ffdBarData;
  
  // Blank out all the Precache variables before reassignment
  PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
  PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
  PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
  PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;
    
  // Make copies of all the precalculations frmo the old runner
  if (old.PreCalcnuPrediction)
    PreCalcnuPrediction = new NuMatrixSpectrum(*old.PreCalcnuPrediction);
  if (old.PreCalcbarPrediction)
    PreCalcbarPrediction = new NuMatrixSpectrum(*old.PreCalcbarPrediction);
  if (old.PreCalcpotentialNuTaus)
    PreCalcpotentialNuTaus = new NuMatrixSpectrum(*old.PreCalcpotentialNuTaus);
  if (old.PreCalcpotentialTauBars)
    PreCalcpotentialTauBars = new NuMatrixSpectrum(*old.PreCalcpotentialTauBars);
  
  if (old.PreCalcnuNCBackground)
    PreCalcnuNCBackground = new NuMatrixSpectrum(*old.PreCalcnuNCBackground);
  if (old.PreCalcbarNCBackground)
    PreCalcbarNCBackground = new NuMatrixSpectrum(*old.PreCalcbarNCBackground);

  if (old.PreCalcwrongSignNuMus)
    PreCalcwrongSignNuMus = new NuMatrixSpectrum(*old.PreCalcwrongSignNuMus);
  if (old.PreCalcwrongSignNuBars)
    PreCalcwrongSignNuBars = new NuMatrixSpectrum(*old.PreCalcwrongSignNuBars);
  
  fPredictNeutrinos = old.fPredictNeutrinos;
  fPredictAntiNeutrinos = old.fPredictAntiNeutrinos;
  
  // Set the wrong sign to false however - this copy may want different
  // oscillation parameters!
  fDoneWsNeutrinos = fDoneWsAntineutrinos = false;
}

//____________________________________________________________________72

NuMMRunFC::~NuMMRunFC()
{
  ResetCache();
}


//___________________________________________________________________

void NuMMRunFC::ResetCache()
{
  // delete the precalculated MM spectra we previously created
  if (PreCalcnuPrediction)    delete PreCalcnuPrediction;
  if (PreCalcbarPrediction)   delete PreCalcbarPrediction;
  
  if (PreCalcpotentialNuTaus) delete PreCalcpotentialNuTaus; 
  if (PreCalcpotentialTauBars)delete PreCalcpotentialTauBars; 
  
  if (PreCalcnuNCBackground)  delete PreCalcnuNCBackground; 
  if (PreCalcbarNCBackground) delete PreCalcbarNCBackground; 
  
  if (PreCalcwrongSignNuMus)  delete PreCalcwrongSignNuMus; 
  if (PreCalcwrongSignNuBars) delete PreCalcwrongSignNuBars; 
  
  PreCalcnuPrediction    = PreCalcbarPrediction    = 0;
  PreCalcpotentialNuTaus = PreCalcpotentialTauBars = 0;
  PreCalcnuNCBackground  = PreCalcbarNCBackground  = 0;
  PreCalcwrongSignNuMus  = PreCalcwrongSignNuBars  = 0;
  
  fCached = false;
}

//____________________________________________________________________72

const std::vector<TH1D>
NuMMRunFC::WriteFDPredHistos(const NuMMParameters& pars) const
{    
  //Set up a vector to push the histograms into.
  vector<TH1D> vHistos;

  //Put the ND data in the vector:
  TH1D hNDNuData(*(fndNuData->Spectrum()));
  hNDNuData.SetName("ndDataNQ");
  hNDNuData.SetTitle("ndDataNQ");
  vHistos.push_back(*(new TH1D(hNDNuData)));
  TH1D hNDBarData(*(fndBarData->Spectrum()));
  hNDBarData.SetName("ndDataPQ");
  hNDBarData.SetTitle("ndDataPQ");
  vHistos.push_back(*(new TH1D(hNDBarData)));

  if (!fQuietMode) {
  cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
     << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
     << endl;
  }

  // Make copies of the prepared extrapolation variables
  NuMatrixSpectrum nuPrediction(*PreCalcnuPrediction);
  NuMatrixSpectrum barPrediction(*PreCalcbarPrediction);
  
  NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTaus);
  NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBars);

  NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
  NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);
  
  NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMus);
  NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBars);
  
  cout << "NC Integral: " << barNCBackground.Spectrum()->Integral() << endl;



// 
//   //Add in backgrounds
//   nuPrediction.Add(wrongSignNuBars);
//   barPrediction.Add(wrongSignNuMus);
//   nuPrediction.Add(nuNCBackground);
//   nuPrediction.Add(potentialNuTaus);
//   barPrediction.Add(barNCBackground);
//   barPrediction.Add(potentialTauBars);
  
  // Processes for neutrinos
  
  if (fPredictNeutrinos)
  {  
    // Oscillate the taus
    if (fHelper->DoTaus()) {
      potentialNuTaus.InverseOscillate(pars.Dm2(),pars.Sn2());
    
      // Write out the true tau spectra
      TH1D hFDTausNuMu(*(potentialNuTaus.Spectrum()));
      hFDTausNuMu.SetName("fdTausNuMu");
      vHistos.push_back(hFDTausNuMu);
      
      potentialNuTaus.TrueToReco(fHelper->FDNuTauRecoVsTrue());
    }

    // Antineutrinos in the neutrino spectrum
    if (!fDoneWsAntineutrinos) wrongSignNuBars.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());

    // Now convert wrong sign to reco
    if (!fDoneWsAntineutrinos) wrongSignNuBars.TrueToReco(fHelper->FDWrongSignBarRecoVsTrue());
    
    // Write out the true wrong sign components
    TH1D hFDWrongSignNQ(*(wrongSignNuBars.Spectrum()));
    hFDWrongSignNQ.SetName("fdWrongSignNQ");
    vHistos.push_back(hFDWrongSignNQ);
    
    // Oscillation and efficiency for the Neutrino spectrum
    nuPrediction.Oscillate(pars.Dm2(),pars.Sn2());
//    nuPrediction.Multiply(fHelper->FDNuEfficiency());

    // Write out the true spectra for the prediction components
    // before we add in backgrounds
    TH1D hFDNoBackTrueNQ(*(nuPrediction.Spectrum()));
    hFDNoBackTrueNQ.SetName("fdBasePredictionNuMus");
    vHistos.push_back(hFDNoBackTrueNQ);

    // True to Reco
    nuPrediction.TrueToReco(fHelper->FDNuRecoVsTrue());

    // Dump out the reco spectrum before adding in the backgrounds
    TH1D hFDNoBackNQ(*(nuPrediction.Spectrum()));
    hFDNoBackNQ.SetName("fdBasePredictionNQ");
    vHistos.push_back(hFDNoBackNQ);

    //Add in backgrounds
    nuPrediction.Add(wrongSignNuBars);
    nuPrediction.Add(nuNCBackground);
    if (fHelper->DoTaus()) nuPrediction.Add(potentialNuTaus);
    
    //Put the other componants of the predictions in the vector
    TH1D hFDPredictionNQ(*(nuPrediction.Spectrum()));
    hFDPredictionNQ.SetName("fdPredictionNQ");
    vHistos.push_back(hFDPredictionNQ);

    TH1D hFDNCNQ(*(nuNCBackground.Spectrum()));
    hFDNCNQ.SetName("fdNCNQ");
    vHistos.push_back(hFDNCNQ);

    if (fHelper->DoTaus()) {
      TH1D hFDTausNQ(*(potentialNuTaus.Spectrum()));
      hFDTausNQ.SetName("fdTausNQ");
      vHistos.push_back(hFDTausNQ);
    }
    
    TH1D hFDDataNQ(*(ffdNuData->Spectrum()));
    hFDDataNQ.SetName("fdDataNQ");
    vHistos.push_back(hFDDataNQ);
    
  } else {
      // Scale the histogram to zero, so as not to confuse people
      nuPrediction.Spectrum()->Scale(0);
  }
  
  // Processes for Antineutrinos
  
  if (fPredictAntiNeutrinos)
  {
    if (fHelper->DoTaus()) {
      // Oscillate the taubars
      potentialTauBars.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
      // Write out the True Tau Spectrum
      TH1D hFDTausNuBar(*(potentialTauBars.Spectrum()));
      hFDTausNuBar.SetName("fdTausNuBar");
      vHistos.push_back(hFDTausNuBar);
      // Convert taus to reco energy
      potentialTauBars.TrueToReco(fHelper->FDTauBarRecoVsTrue());
    }
    // Oscillate the wrong sign component
    // Neutrinos for the antineutrino spectrum
    if (!fDoneWsNeutrinos) wrongSignNuMus.Oscillate(pars.Dm2(),pars.Sn2());

    // Now, convert wrong sign to reco energy
    if (!fDoneWsNeutrinos) wrongSignNuMus.TrueToReco(fHelper->FDWrongSignNuRecoVsTrue());
    
    // write out the true wrong sign spectra
    TH1D hFDWrongSignPQ(*(wrongSignNuMus.Spectrum()));
    hFDWrongSignPQ.SetName("fdWrongSignPQ");
    vHistos.push_back(hFDWrongSignPQ);
    
    //Oscillation and efficiency for the antineutrino spectrum
    barPrediction.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());

    // Write out the true antineutrino spectrum, without backgrounds
    TH1D hFDNoBackTruePQ(*(barPrediction.Spectrum()));
    hFDNoBackTruePQ.SetName("fdBasePredictionNuBars");
    vHistos.push_back(hFDNoBackTruePQ);

    //True to reco
    barPrediction.TrueToReco(fHelper->FDBarRecoVsTrue());

    // Antineutrino reco spectrum, before backgrounds
    TH1D hFDNoBackPQ(*(barPrediction.Spectrum()));
    hFDNoBackPQ.SetName("fdBasePredictionPQ");
    vHistos.push_back(hFDNoBackPQ);

    //Add in backgrounds
    barPrediction.Add(wrongSignNuMus);
    barPrediction.Add(barNCBackground);
    if (fHelper->DoTaus()) barPrediction.Add(potentialTauBars);
    
    TH1D hFDPredictionPQ(*(barPrediction.Spectrum()));
    hFDPredictionPQ.SetName("fdPredictionPQ");
    vHistos.push_back(hFDPredictionPQ);

    TH1D hFDNCPQ(*(barNCBackground.Spectrum()));
    hFDNCPQ.SetName("fdNCPQ");
    vHistos.push_back(hFDNCPQ);

    if (fHelper->DoTaus()) {
      TH1D hFDTausPQ(*(potentialTauBars.Spectrum()));
      hFDTausPQ.SetName("fdTausPQ");
      vHistos.push_back(hFDTausPQ);
    }

    TH1D hFDDataPQ(*(ffdBarData->Spectrum()));
    hFDDataPQ.SetName("fdDataPQ");
    vHistos.push_back(hFDDataPQ);
    
  } else {
      // Scale the histogram to zero, so as not to confuse people
      barPrediction.Spectrum()->Scale(0);
  }

  return vHistos;
}


//___________________________________________________________________

void NuMMRunFC::CacheExtrapolation( void )
{
    // Only cache once.
    if (fCached) return;
  
    ResetCache();
  
    // Disable this output for now - this process should be fast anyway!
    // if (!fQuietMode)
    //     cout << "Pre-calculating static extrapolation variables..." << endl;
    
    // Caching variables, to reduce the workload of each step of the fitter
    // Initialise these the way they would be before
    PreCalcnuPrediction = new NuMatrixSpectrum(*fndNuData);
    PreCalcbarPrediction = new NuMatrixSpectrum(*fndBarData);
    
    // Be lazy and just alias the pointers to references so can keep
    // the code the same
    NuMatrixSpectrum &nuPrediction = *PreCalcnuPrediction;
    NuMatrixSpectrum &barPrediction = *PreCalcbarPrediction;
    
    // Now we have prepared the extrapolation variables, calculate them
    // Get the neutrinos to the FD
    nuPrediction.Multiply(fHelper->NDNuPurity());
    nuPrediction.RecoToTrue(fHelper->NDNuRecoVsTrue());
    nuPrediction.Divide(fHelper->NDNuEfficiency());
    nuPrediction.Divide(fHelper->NuXSecGraph());
    nuPrediction.Divide(fndNuData->PoT());
    nuPrediction.Divide(fNDFidMass);
    nuPrediction.ExtrapolateNDToFD(fHelper->NuBeamMatrix());
    nuPrediction.Multiply(fFDFidMass);
    nuPrediction.Multiply(ffdNuData->PoT());
    // Reset the spectrum to the new PoT
    nuPrediction.ResetPoT(ffdNuData->PoT());

    //Get the antineutrinos to the FD
    barPrediction.Multiply(fHelper->NDBarPurity());
    barPrediction.RecoToTrue(fHelper->NDBarRecoVsTrue());
    barPrediction.Divide(fHelper->NDBarEfficiency());
    barPrediction.Divide(fHelper->BarXSecGraph());
    barPrediction.Divide(fndBarData->PoT());
    barPrediction.Divide(fNDFidMass);
    barPrediction.ExtrapolateNDToFD(fHelper->BarBeamMatrix());
    barPrediction.Multiply(fFDFidMass);
    barPrediction.Multiply(ffdBarData->PoT());
    // Reset the spectrum to the new PoT
    barPrediction.ResetPoT(ffdBarData->PoT());

    // Create and link to the Tau components, before we include the
    // cross section in the predictions (tau Xsec is different)
    PreCalcpotentialNuTaus = new NuMatrixSpectrum(nuPrediction);
    PreCalcpotentialTauBars = new NuMatrixSpectrum(barPrediction);
    NuMatrixSpectrum &potentialNuTaus = *PreCalcpotentialNuTaus;
    NuMatrixSpectrum &potentialTauBars = *PreCalcpotentialTauBars;

    if (fHelper->DoTaus()) {
      //Get the taus (ignoring any wrong-sign taubars)
      potentialNuTaus.Multiply(fHelper->XSecGraphNuTaus());
      potentialNuTaus.Multiply(fHelper->FDNuTauEfficiency());
      //Get the taubars (ignoring any wrong-sign taus)
      potentialTauBars.Multiply(fHelper->XSecGraphTauBars());
      potentialTauBars.Multiply(fHelper->FDTauBarEfficiency());
    }
    else {
      PreCalcpotentialNuTaus->Multiply(0.0);
      PreCalcpotentialTauBars->Multiply(0.0);
    }
  
    //Cross sections
    nuPrediction.Multiply(fHelper->NuXSecGraph());
    barPrediction.Multiply(fHelper->BarXSecGraph());

    // Now build the NC Spectrum, after the cross section calculation!
    PreCalcnuNCBackground = new NuMatrixSpectrum(*PreCalcnuPrediction);;
    PreCalcbarNCBackground = new NuMatrixSpectrum(*PreCalcbarPrediction);
    NuMatrixSpectrum &nuNCBackground = *PreCalcnuNCBackground;
    NuMatrixSpectrum &barNCBackground = *PreCalcbarNCBackground;  

    //Get the neutrino NC background
    nuNCBackground.Multiply(fHelper->FDNuEfficiency());
    nuNCBackground.TrueToReco(fHelper->FDNuRecoVsTrue());
    nuNCBackground.Divide(fHelper->FDNuPurity());
    nuNCBackground.Multiply(fHelper->FDNuNCContamination());

    //Get the antineutrino NC background
    barNCBackground.Multiply(fHelper->FDBarEfficiency());
    barNCBackground.TrueToReco(fHelper->FDBarRecoVsTrue());
    barNCBackground.Divide(fHelper->FDBarPurity());
    barNCBackground.Multiply(fHelper->FDBarNCContamination());
    
    // Now build the wrong sign spectrum
    PreCalcwrongSignNuMus = new NuMatrixSpectrum(*PreCalcnuPrediction);
    PreCalcwrongSignNuBars = new NuMatrixSpectrum(*PreCalcbarPrediction);
    NuMatrixSpectrum &wrongSignNuMus = *PreCalcwrongSignNuMus;
    NuMatrixSpectrum &wrongSignNuBars = *PreCalcwrongSignNuBars;
    
    // Do the Neutrino background for the antineutrino spectrum
    wrongSignNuMus.Multiply(fHelper->FDWrongSignNuEfficiency());
    
    // Do the Antineutrino background for the neutrino spectrum
    wrongSignNuBars.Multiply(fHelper->FDWrongSignBarEfficiency());
    
    // Apply the efficiency to the predictions
    barPrediction.Multiply(fHelper->FDBarEfficiency());
    nuPrediction.Multiply(fHelper->FDNuEfficiency());
    
  
    nuPrediction.SetName("NuMus");
    barPrediction.SetName("NuMuBars");
    
    if (fHelper->DoTaus()) {
      potentialNuTaus.SetName("NuTaus");
      potentialTauBars.SetName("NuTauBars");
    }
    
    wrongSignNuMus.SetName("WrongSign_NuMus");
    wrongSignNuBars.SetName("WrongSign_NuMuBars");
    
    nuNCBackground.SetName("NC_NuMus");
    barNCBackground.SetName("NC_NuMuBars");
  
    fCached = true;
  
    // If we have precalculated the wrong sign oscillations,
    // then we need to recalculate them next time we are given
    // the oscilaltion parameters. Reset the flags used to
    // control this  
    fDoneWsNeutrinos = fDoneWsAntineutrinos = false;
}

//____________________________________________________________________

void NuMMRunFC::PreCalcWrongSign(const NuMMParameters& pars)
{
    // Precalculate the wrong sign, if we have disabled predictions
    // for that class, and it hasn't been done already. This should only
    // be done once, unless of course the user does something strange like
    // turning on neutrino prediction in the middle of an extrapolation.
    if(!fDoneWsNeutrinos) {
        if (!fPredictNeutrinos)
        {
            if (!fQuietMode) cout << "Neutrino prediction off: Precalculating wrong sign oscillation." << endl;
            PreCalcwrongSignNuMus->Oscillate(pars.Dm2(),pars.Sn2());
            PreCalcwrongSignNuMus->TrueToReco(fHelper->FDWrongSignNuRecoVsTrue());
            fDoneWsNeutrinos = true;
        }
    }
    if(!fPredictAntiNeutrinos) {
        if (!fDoneWsAntineutrinos)
        {
            if (!fQuietMode) cout << "Antineutrino prediction off: Precalculating wrong sign oscillation." << endl;
            PreCalcwrongSignNuBars->Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
            PreCalcwrongSignNuBars->TrueToReco(fHelper->FDWrongSignBarRecoVsTrue());
            fDoneWsAntineutrinos = true;
        }
    }
}

//____________________________________________________________________

const pair<NuMatrixSpectrum,NuMatrixSpectrum>
NuMMRunFC::MakeFDPred(const NuMMParameters& pars)
{
  // Precache extrap, but only if needed
  CacheExtrapolation(); 
  // Do the wrong sign precalcuation, if it hasn't been done
  PreCalcWrongSign(pars);
  
  // Output the oscillation settings!
  if (!fQuietMode) {
      //pars.PrintStatus();
    cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
    << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
    << endl;
  }
  
  // For all neutrinos

  // Make copies of the prepared extrapolation variables, so that
  // we can alter them locally with the oscillation parameters
  NuMatrixSpectrum nuPrediction(*PreCalcnuPrediction);
  NuMatrixSpectrum barPrediction(*PreCalcbarPrediction);
  
  NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTaus);
  NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBars);

  NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
  NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);
  
  NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMus);
  NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBars);

  // Processes for neutrinos
  
  if (fPredictNeutrinos)
  {  
    if (fHelper->DoTaus()) {
      // Oscillate the taus
      potentialNuTaus.InverseOscillate(pars.Dm2(),pars.Sn2());
      potentialNuTaus.TrueToReco(fHelper->FDNuTauRecoVsTrue());
    }

    // If we have pre-oscillated the wrong sign, do not do so again
    if (!fDoneWsAntineutrinos) {
        // Antineutrinos in the neutrino spectrum
        wrongSignNuBars.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
        wrongSignNuBars.TrueToReco(fHelper->FDWrongSignBarRecoVsTrue());
    }

    // Oscillate the Neutrino spectrum!
    nuPrediction.Oscillate(pars.Dm2(),pars.Sn2());

    // True to Reco
    nuPrediction.TrueToReco(fHelper->FDNuRecoVsTrue());

    //Add in backgrounds
    nuPrediction.Add(wrongSignNuBars);
    nuPrediction.Add(nuNCBackground);
    if (fHelper->DoTaus()) nuPrediction.Add(potentialNuTaus);
  } else {
      // Scale the histogram to zero, so as not to confuse people
      // if they look at the spectrum with neutrinos turned off
      nuPrediction.Spectrum()->Scale(0);
  }
  
  // Processes for Antineutrinos
  
  if (fPredictAntiNeutrinos)
  {
      if (fHelper->DoTaus()) {
        // Oscillate the taubars
        potentialTauBars.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
        potentialTauBars.TrueToReco(fHelper->FDTauBarRecoVsTrue());
      }

      // Only oscillate the wrong sign if it hasn't been done already
      if (!fDoneWsNeutrinos) {
          // Oscillate the wrong sign component
          // Neutrinos for the antineutrino spectrum
          wrongSignNuMus.Oscillate(pars.Dm2(),pars.Sn2());
          wrongSignNuMus.TrueToReco(fHelper->FDWrongSignNuRecoVsTrue());
      }
      // Oscillate the antineutrino spectrum!
      barPrediction.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());

      // True to reco
     barPrediction.TrueToReco(fHelper->FDBarRecoVsTrue());
      
      // Add in backgrounds
      barPrediction.Add(wrongSignNuMus);
      barPrediction.Add(barNCBackground);
      if (fHelper->DoTaus()) barPrediction.Add(potentialTauBars);
  } else {
      // Scale the histogram to zero, so as not to confuse people
      barPrediction.Spectrum()->Scale(0);
  }
  
  // Return the two predictions, as a pair
  pair<NuMatrixSpectrum,NuMatrixSpectrum>
    predictions(nuPrediction,barPrediction);
  return predictions;
}


//____________________________________________________________________
// Be careful -- DO NOT precalculate the wrong signs for this function

const NuMatrixSpectrum NuMMRunFC::TrueComponents(const NuMMParameters& pars, Sample_t s)
{
  // Precache extrap, but only if needed
  CacheExtrapolation(); 

  // Output the oscillation settings!
  if (!fQuietMode) {
    cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
    << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
    << endl;
  }
  
  if (fDoneWsNeutrinos || fDoneWsAntineutrinos) {
    MSG("NuMMRunFC",Msg::kError) << "Already True-To-Recoed the wrong signs.  Cannot precalculate the wrong-signs for event-by-event.  Bailing." << endl;
    MSG("NuMMRunFC",Msg::kFatal) << "Already True-To-Recoed the wrong signs.  Cannot precalculate the wrong-signs for event-by-event.  Bailing." << endl;
  } 
  
  // No appeared spectrum for oscillations
  if (s == kAppearedPQ || s == kAppearedNQ) {
    NuMatrixSpectrum blank(*PreCalcnuPrediction);
    blank.SetName("Blank");
    blank.Spectrum()->Scale(0);
    return blank;    
  }
  
  // Request NQ spectrum, but no NuBar Predictions OR
  // Request PQ spectrum, but no NuMu Predictions
  if ( (s < 0 && !fPredictNeutrinos) || 
      (s > 0 && !fPredictAntiNeutrinos) ) {
    NuMatrixSpectrum blank(*PreCalcnuPrediction);
    blank.SetName("Blank");
    blank.Spectrum()->Scale(0);
    return blank;
  }
  
  if (s == kSignalPQ) {
    NuMatrixSpectrum barPrediction(*PreCalcbarPrediction);
    barPrediction.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
    return barPrediction;
  }
  else if (s == kSignalNQ) {
    NuMatrixSpectrum nuPrediction(*PreCalcnuPrediction);  
    nuPrediction.Oscillate(pars.Dm2(),pars.Sn2());
    return nuPrediction;
  }
  else if (s == kWrongSignPQ) {
    NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMus);
    wrongSignNuMus.Oscillate(pars.Dm2(),pars.Sn2());
    return wrongSignNuMus;
  }
  else if (s == kWrongSignNQ) {
    NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBars);
    wrongSignNuBars.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
    return wrongSignNuBars;
  }  
  else if (s == kNCPQ) {
    NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);
    return barNCBackground;
  }
  else if (s == kNCNQ) {
    NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
    return nuNCBackground;
  }
  else if (s == kTauPQ) {
    NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBars);
    potentialTauBars.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
    return potentialTauBars;
  }
  else if (s == kTauNQ) {
    NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTaus);
    potentialNuTaus.InverseOscillate(pars.Dm2(),pars.Sn2());
    return potentialNuTaus;
  }
  else {
    MSG("NuMMRunTransitions",Msg::kWarning) << "Requested unknown prediction #" << s << ", returning blank." << endl;
    NuMatrixSpectrum blank(*PreCalcnuPrediction);
    blank.SetName("Blank");
    blank.Spectrum()->Scale(0);
    return blank;
  }
}

//____________________________________________________________________

const Double_t NuMMRunFC
::ComparePredWithData(const NuMMParameters& pars)
{
  // Generate a prediction
  const pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions =
    this->MakeFDPred(pars);
  Double_t like = 0;
  // Now, only count a likelihood for components we want - this is just
  // a constant offset, but should help speed things up a little.
  if (fPredictNeutrinos)
  {
    like += this->StatsLikelihood(predictions.first.Spectrum(),
        ffdNuData->Spectrum());
  }
  if (fPredictAntiNeutrinos) {
    like += this->StatsLikelihood(predictions.second.Spectrum(),
        ffdBarData->Spectrum());
  }
  return like;
}

---