NuFCExperimentFactory.cxx

Go to the documentation of this file.

#include <string>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <stdexcept>
#include <map>

using namespace std;

#include "TROOT.h"
#include "TH1.h"
#include "TH1D.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TFile.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuFCExperimentFactory.h"
#include "NtupleUtils/NuMMParameters.h"
#include "NtupleUtils/NuInputEvents.h"
#include "NtupleUtils/NuLibrary.h"
#include "NtupleUtils/NuFCEvent.h"
#include "NtupleUtils/NuMMRunNuBar.h"


CVSID("$Id: NuFCExperimentFactory.cxx,v 1.26 2010/01/26 11:56:23 bckhouse Exp $");

using namespace Sample;

NuFCExperimentFactory::NuFCExperimentFactory(TRandom3 *rgen, TString nd_path, TString fd_path, TString tau_path)
  : IntND(0), ndPoT(0), 
    fPQHist(), fNQHist(), fPQHistND(), fNQHistND(), 
    fRandy(rgen), 
    counter(0), 
    fixedSyst(false), gausSyst(false), dp_corr(false), dp_old(false)
{
  // Check if our random number generator is correctly setup
  if (!fRandy) {
    // If not, just use gRandom
    fRandy = gRandom;
  }
  noND = false;
  
  norm_size    = 0.04;

  curv_size    = 0.02;
  relcurv_size = 0.04;
  range_size   = 0.02;
  shower_size  = 0.10;
  relshow_size = 0.033;
  
  dp_size      = 0.80;
  bg_size      = 0.50;
  
  xsec_size    = 1;
  skzp_size    = 1;

  // Get the binning scheme set up
  const NuUtilities cuts;
  NuBinningScheme::NuBinningScheme_t binningScheme =
  static_cast<NuBinningScheme::NuBinningScheme_t>(4);
  vReco = cuts.RecoBins(binningScheme);
  vTrue = cuts.TrueBins(binningScheme);
  
  numRecoBins = vReco.size() - 1;
  numTrueBins = vTrue.size() - 1;
  //recoBinsArray = &(vReco[0]);
  //trueBinsArray = &(vTrue[0]);

  Distro(kSignalPQ, "signalSpect", "Signal Distro");
  Distro(kWrongSignPQ, "wsSpect", "Wrong Sign Distro");
  Distro(kNCPQ, "ncSpect", "NC Distro");
  Distro(kTauPQ, "tauSpect", "Tau Distro");
  Distro(kAppearedPQ, "appearedSpect", "Appeared Distro");
  
  // Set up the FD MC chain
  FillMCEvents(nd_path);
  FillMCEvents(fd_path);
  FillMCEvents(tau_path);
}


void NuFCExperimentFactory::GenerateNewExperiment(NuMMRunNuBar * mmRunSource, const NuMMParameters& mmPars, bool reroll)
{
  MAXMSG("NuFCExperimentFactory",Msg::kInfo,1) << "Entering GenerateNewExperiment" << endl;
  
  TString name = "fdExp";
  name += (counter++);
  TH1D *fdExp = new TH1D(name,"Single FC Experiment", numRecoBins, &(vReco[0]));

  NuMatrixSpectrum fdBar;
  mmRunSource->QuietModeOn();

  if (reroll) RollSystematics();
  GenerateND(mmRunSource);
  
  for (int i = 1; i <= 5; i++) {
    Sample_t s = (Sample_t)i;
    
    fdBar = mmRunSource->TrueComponents(mmPars, s);
    
    // Apply Normalization to the Prediction
    fdBar.Multiply(GetNormShift());
    fdBar.Multiply(GetBGShift(s));
    
    // Get the average events by integrating the PDF
    Double_t average_events = fdBar.Spectrum()->Integral();
    counts[Int(s)] = 0;
    if (average_events == 0) continue;
    
    // Get the number of events to generate
    Int_t nevents = fRandy->Poisson(average_events);

    // Create sampling distribution
    MakeSampler(s, fdBar);
    
    double numdraws = 0;

    // Now start sampling
    int ev;
    double rand, reject, E;
    bool seeking;
    
    while (nevents > 0) {
      // Get an event -- takes care of rejection sampling
      seeking = true;
      while (seeking) {
        ev = fRandy->Integer(N(s));
        ++numdraws;
        
        rand = fRandy->Uniform();
        reject = GetReject(s, ev);       
        
        if (s != kAppearedPQ) {
          ScaleRejection(reject, nevents, rand, GetDPShift(events[Int(s)][ev]));
        }
        if (s != kWrongSignPQ) {
          ScaleRejection(reject, nevents, rand, GetXSecShift(events[Int(s)][ev]));
        }
        ScaleRejection(reject, nevents, rand, GetSKZPShift(events[Int(s)][ev]));

        
        seeking = rand > reject;
      }

      // Get the new reco energy -- apply the systematic
      E = GetTrkEnergy(events[Int(s)][ev]) + GetShowerEnergy(events[Int(s)][ev]);

      // Fill our final reco spectrum
      fdExp->Fill(E);
      ++counts[Int(s)];
      
      --nevents;
    }
  }
  
  fNQHist.ResetPoT(fdBar.PoT());
  fPQHist.ResetPoT(fdBar.PoT());
  fPQHist.ResetSpectrum(*fdExp);

  delete fdExp;
  
  MAXMSG("NuFCExperimentFactory",Msg::kInfo,1) << "Leaving GenerateNewExperiment" << endl;
}


void NuFCExperimentFactory::ScaleRejection(double &reject, int &nevents, double rand, double shift)
{ 
  double reject_new = reject*shift;
  if (rand > reject && rand < reject_new) nevents++;
  if (rand < reject && rand > reject_new) nevents--;
  
  reject = reject_new;
}


void NuFCExperimentFactory::FillMCEvents(TString mc_path)
{
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Entering FillMCEvents(\"" << mc_path << "\")" << endl;
  TFile *fin = new TFile(mc_path,"READ");

  // Test the file opened, and fail otherwise
  if (fin->IsZombie()) {
    MSG("NuFCExperimentFactory",Msg::kError) << "Failed to read in events file " << mc_path << endl;
    throw runtime_error("Failed to read MC events file");
  }
  
  // Determine which Detector this file is for
  TH1D *hDetector = (TH1D*)gROOT->FindObject("hDetector");
  int det = (int)hDetector->GetMean();
 
  // Get the ND Pot
  if (det == 1) {
    TH1D *hTotalPot = (TH1D*)gROOT->FindObject("hTotalPot");
    ndPoT += hTotalPot->Integral();
  }
  
  // Attempt to load the known tree types
  TTree *t1 = (TTree*)fin->Get("FDFilter");
  TTree *t2 = (TTree*)fin->Get("TAUFilter");
  TTree *t3 = (TTree*)fin->Get("NDFilter");
  TTree *t4 = (TTree*)fin->Get("FCTree");
  TTree *t5 = (TTree*)fin->Get("s");

  // If the tree exists in the file, load it
  if (t1)      FillFromFilter(t1, det);
  else if (t2) FillFromFilter(t2, det);
  else if (t3) FillFromFilter(t3, det);
  else if (t4) FillFromFCTree(t4, det);
  else if (t5) FillFromDST(mc_path, det);
  else {
    MSG("NuFCExperimentFactory",Msg::kError) << "No recognized tree in file: " << mc_path << endl;      
  }
  
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Leaving FillMCEvents" << endl;  
}  


void NuFCExperimentFactory::FillFromFilter(TTree *tin, int det)
{
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Entering FillFromFilter" << endl;
  
  Float_t energy;       tin->SetBranchAddress("energy",       &energy);
  Float_t trkEn;        tin->SetBranchAddress("trkEn",        &trkEn);
  Float_t shwEn;        tin->SetBranchAddress("shwEn",        &shwEn);
  Float_t neuEnMC;      tin->SetBranchAddress("neuEnMC",      &neuEnMC);  
  Int_t   iaction;      tin->SetBranchAddress("iaction",      &iaction);    
  Int_t   inu;          tin->SetBranchAddress("inu",          &inu);
  Float_t charge;       tin->SetBranchAddress("charge",       &charge);
  Float_t ppvz;         tin->SetBranchAddress("ppvz",         &ppvz); 
  Int_t   containedTrk; tin->SetBranchAddress("containedTrk", &containedTrk); 
  Float_t fluxErr;      tin->SetBranchAddress("fluxErr",      &fluxErr); 
  
  float vals[10];
  for (Int_t i=0; i < tin->GetEntries(); ++i) {
    tin->GetEntry(i);
    this->PrintLoopProgress(i,tin->GetEntries());
   
    vals[0] = energy;
    vals[1] = trkEn;
    vals[2] = shwEn;
    vals[3] = neuEnMC;
    vals[4] = (float)iaction;
    vals[5] = (float)inu;
    vals[6] = charge;
    vals[7] = ppvz;    
    vals[8] = (float)containedTrk;
    vals[9] = fluxErr;
    
    if (inu == 16 || inu == -16) {
      // Ignore taus from Nue's, NC's
//      if (nu->inunoosc == 12 || nu->inunoosc==-12 || nu->iaction == 0) continue;
      if (iaction == 0) continue;
      StoreEvent(kTauPQ, vals, det);
    }
    else if (iaction == 0) StoreEvent(kNCPQ, vals, det);
    else if (inu == -14)   {
      StoreEvent(kAppearedPQ, vals, det);
      StoreEvent(kSignalPQ, vals, det);
    }
    else StoreEvent(kWrongSignPQ, vals, det);
  }
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Leaving FillFromFilter" << endl;  
}
 

void NuFCExperimentFactory::FillFromFCTree(TTree *tin, int det)
{
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Entering FillFromFCTree" << endl;
  
  NuFCEvent *nu = new NuFCEvent();
  tin->SetBranchAddress("NuFCEvent", &nu);

  for (Int_t i=0; i < tin->GetEntries(); ++i) {
    tin->GetEntry(i);
    this->PrintLoopProgress(i,tin->GetEntries());

    if (nu->inu == 16 || nu->inu == -16) {
      // Ignore taus from Nue's, NC's
      if (nu->inunoosc == 12 || nu->inunoosc==-12 || nu->iaction == 0) continue;
      StoreEvent(kTauPQ, *nu, det);
    }
    else if (nu->iaction == 0) StoreEvent(kNCPQ, *nu, det);
    else if (nu->inu == -14)   {
      StoreEvent(kAppearedPQ, *nu, det);
      StoreEvent(kSignalPQ, *nu, det);
    }
    else StoreEvent(kWrongSignPQ, *nu, det);
  }
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Leaving FillFromFCTree" << endl;  
}


void NuFCExperimentFactory::FillFromDST(TString mc_path, int det)
{
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Entering FillFromDST(" << mc_path << ")" << endl;
  
  // TChain holding all the FD MC
  NuInputEvents* mc_input = new NuInputEvents();
  mc_input->InputFileName(mc_path.Data());
  mc_input->InitialiseChains();
  mc_input->InitialiseNuEventBranch();  
  mc_input->ResetNuEventLoopPositionToStart();
  
  for (Int_t i=0; i < mc_input->GetEntriesNuEvent(); ++i) {
    NuEvent &nu=const_cast<NuEvent&>(mc_input->GetNextNuEvent(Msg::kInfo));
    
    //this->PrintLoopProgress(i,mc_input->GetEntriesNuEvent());
    
    if (!IsGoodStdCuts(nu)) continue;
    
    if (nu.inu == 16 || nu.inu == -16) {
      // Ignore taus from Nue's, NC's
      if (nu.inunoosc == 12 || nu.inunoosc==-12 || nu.iaction == 0) continue;
      StoreEvent(kTauPQ, nu, det);
    }
    else if (nu.iaction == 0) StoreEvent(kNCPQ, nu, det);
    else if (nu.inu == -14)   {
      StoreEvent(kAppearedPQ, nu, det);
      StoreEvent(kSignalPQ, nu, det);
    }
    else StoreEvent(kWrongSignPQ, nu, det);
  }
  
  MSG("NuFCExperimentFactory",Msg::kInfo) << "Leaving FillFromDST" << endl;  
}


Bool_t NuFCExperimentFactory::IsGoodStdCuts(const NuEvent& nu) const
{
  // Only keep reconstructed NuBars
  if (nu.charge == -1) return false;
    
  //get an instance of the code library
  NuLibrary& lib=NuLibrary::Instance();
  
  //cut on the sntp good beam and that coil is on
  if (!lib.cuts.IsGoodBeamDetPOTCountingStage(nu)) return false;
  //lib.cnt.goodBeamDetPOTCountingStage++;
  
  //ensure good number of tracks in the event
  if (!lib.cuts.IsGoodNumberOfTracks(nu)) return false;
  lib.cnt.evtWithTrkCounter++;
  
  //check if the trk is in the fiducial volume
  if (!lib.cuts.IsInFidVolTrk(nu)) return false;
  lib.cnt.trkInFidVolCounter++;
  
  //cut on LI
  if (lib.cuts.IsLI(nu)) return false;
  lib.cnt.evtNotIsLI++;
  
  //cut on the data quality
  if (!lib.cuts.IsGoodDataQuality(nu)) return false;
  lib.cnt.goodDataQualityCounter++;
  
  //cut on the spill time
  if (!lib.cuts.IsGoodTimeToNearestSpill(nu)) return false;
  lib.cnt.goodTimeToNearestSpillCounter++;
  
  //cut on the beam
  if (!lib.cuts.IsGoodBeam(nu)) return false;
  lib.cnt.goodBeamInfoDBCounter++;
  
  //require a good trk fit
  if (!lib.cuts.IsGoodTrackFitPass(nu)) return false;
  lib.cnt.goodTrkPassCounter++;
  
  //require a forward going neutrino about beam direction
  if (!lib.cuts.IsGoodDirCos(nu)) return false;
  lib.cnt.goodDirectionCosineCounter++;
  
  //cut on the PID
  if (!lib.cuts.IsGoodPID(nu)) return false;
  lib.cnt.goodPIDCounter++;   
  
  //cut on the fractional track momentum and sign error      
  if (!lib.cuts.IsGoodSigmaQP_QP(nu)) return false;
  lib.cnt.goodFitSigQPCounter++;
  
  //cut on the track fit probability      
  if (!lib.cuts.IsGoodFitProb(nu)) return false;
  lib.cnt.goodFitProbCounter++;    
  
  //cut on majority curvature
  if (!lib.cuts.IsGoodMajorityCurvature(nu)) return false;
  
  // Cut on the track length
  if (!lib.cuts.IsGoodTrackLength(nu)) return false;
  
  // Cut on the relative angle
  if (!lib.cuts.IsGoodRelAngle(nu)) return false;
  
  return true;
}


void NuFCExperimentFactory::GenerateND(NuMMRunNuBar * mmRunSource)
{
  if (noND) {
    MAXMSG("NuFCExperimentFactory",Msg::kInfo,1) << "Not Generating the ND" << endl;
    fPQHistND = *(mmRunSource->GetNDBarData());
  }
  else {  
    MAXMSG("NuFCExperimentFactory",Msg::kInfo,1) << "Entering GenerateND" << endl;
    TString name = "ndSpect";
    name += counter;
    TH1D newND(name,"ND Spect", numRecoBins, &(vReco[0]));
    TH1D newDP(name+"_dp","ND Spect, DP", numRecoBins, &(vReco[0]));      
    /*
     unsigned int max = event[IntND].size();
     if (max == 0) {
     MSG("NuFCExperimentFactor",Msg::kError) << "No ND events loaded.  Cannot make an ND spectrum."
     }
     */
    
    double E;
    double rw;
    vector<NuFCEvent*>::const_iterator it;
    for (it = events[IntND].begin(); it != events[IntND].end(); ++it) {
      E = GetTrkEnergy(*it) + GetShowerEnergy(*it);
      rw = (*it)->rw;
      rw *= GetBGShift(*it);
      rw *= GetXSecShift(*it);
      rw *= GetSKZPShift(*it);
      
      if (dp_corr) {
        if (IsDP(*it)) newDP.AddBinContent(EBin(E), rw);
        else           newND.AddBinContent(EBin(E), rw);
      }
      else {
        rw *= GetDPShift(*it);
        newND.AddBinContent(EBin(E), rw);
      }
    }
    
    if (dp_corr) {
      // Set the decay pipe systematic size
      int lowb = 1;
      int highb = newND.GetXaxis()->FindFixBin(13.0);
      double U = newND.Integral(lowb, highb) / ndPoT * 2.92878e+20;
      double D = newDP.Integral(lowb, highb)/ ndPoT * 2.92878e+20;
      double N = 198064.;
      dp_shift = (N - U) / D;
      if (dp_shift < 0) dp_shift = 0;
      newDP.Scale(dp_shift);
      newND.Add(&newDP);
    }
    
    fPQHistND.ResetSpectrum(newND);
    fPQHistND.ResetPoT(ndPoT);
    MAXMSG("NuFCExperimentFactory",Msg::kInfo,1) << "Leaving GenerateND" << endl;
  }

  fNQHistND = *(mmRunSource->GetNDNuData());
}


int NuFCExperimentFactory::EBin(double E) const
{
  if (E < 0) return 0;
  if (E < 0.5) return 1;
  if (E < 20.) return static_cast<int>(E*4.);
  if (E < 30.) return static_cast<int>(E)+60;
  if (E < 50.) return static_cast<int>(E)/2+75;
  return 100;
}



void NuFCExperimentFactory::PrintLoopProgress(Int_t entry,Float_t nEntries) const
{
  Float_t fract=ceil(nEntries/20.);  
  if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
    MSG("NuFCExperimentFactory",Msg::kInfo) 
    <<"Fraction of loop complete: "<<entry 
    <<"/"<<nEntries<<"  ("
    <<(Int_t)(100.*entry/nEntries)<<"%)"<<endl;
  }
}

  


int NuFCExperimentFactory::Int(Sample_t s) const
{
  switch (s) {
      case kSignalPQ:
        return 1;
      case kWrongSignPQ:
        return 2;
      case kNCPQ:
        return 3;
      case kTauPQ:
        return 4;
      case kAppearedPQ:
        return 5;
      default:
        MSG("NuFCExperimentFactor",Msg::kError) << s << " is not a known Sample_t.  Returning -1" << endl;
        return -1;
  }
}


void NuFCExperimentFactory::Distro(Sample_t s, TString name, TString title) 
{
  spect[Int(s)] = new TH1D(name, title, numTrueBins, &(vTrue[0]));
  hSampler[Int(s)] = new TH1D();  // Prevents seg faults
}


TH1D * NuFCExperimentFactory::Distro(Sample_t s) const
{
  return spect[Int(s)];
}


void NuFCExperimentFactory::StoreEvent(Sample_t s, float *vals, int det)
{
  NuFCEvent loc;
  loc.energy       = vals[0];
  loc.trkEn        = vals[1];
  loc.shwEn        = vals[2];
  loc.neuEnMC      = vals[3];
  loc.iaction      = (Int_t)vals[4];
  loc.inu          = (Int_t)vals[5];
  loc.charge       = (Int_t)vals[6];
  loc.ppvz         = vals[7];
  loc.containedTrk = (Int_t)vals[8];
  loc.fluxErr      = vals[9];    
  
  StoreEvent(s, loc, det);
}


void NuFCExperimentFactory::StoreEvent(Sample_t s, NuEvent ev, int det)
{
  NuFCEvent loc(ev);
  StoreEvent(s, loc, det);
}


void NuFCExperimentFactory::StoreEvent(Sample_t s, NuFCEvent ev, int det)
{
  NuFCEvent * loc = new NuFCEvent(ev);
  loc->detector = det;
  
  if (loc->detector == 2) { // Far Det and Taus
    events[Int(s)].push_back(loc);
    
    double E;
    if (s == kNCPQ) E = loc->energy;
    else            E = loc->neuEnMC;
    
    Distro(s)->Fill(E);
  }
  else if (loc->detector == 1) {
    // Prevent double counting the appeared events.
    if (s != kAppearedPQ) events[IntND].push_back(loc);
  }
  else {
    MSG("NuFCExperimentFactory",Msg::kWarning) << "Detector " << det << " not recognized.  Event not stored." << endl;
  }
}


double NuFCExperimentFactory::GetReject(Sample_t s, int i) const
{
  double Esel, reject;
  if (s == kNCPQ) Esel = events[Int(s)][i]->energy;
  else            Esel = events[Int(s)][i]->neuEnMC;
  
  int bin = hSampler[Int(s)]->GetXaxis()->FindFixBin(Esel);
  reject = hSampler[Int(s)]->GetBinContent(bin);
  return reject;
}


void NuFCExperimentFactory::MakeSampler(Sample_t s, NuMatrixSpectrum pred)
{
  if (hSampler[Int(s)]) delete hSampler[Int(s)];
  hSampler[Int(s)] = new TH1D(*(pred.Spectrum()));
  hSampler[Int(s)]->Divide(Distro(s));
  hSampler[Int(s)]->Scale(1./hSampler[Int(s)]->GetMaximum());
}


void NuFCExperimentFactory::DebugMode(bool deb)
{
  fixedSyst = deb;
  
  if (fixedSyst) {
    MSG("NuFCExperimentFactory",Msg::kInfo) << "Now in debugging mode.  Systematics are not random -- they are fixed at the input values." << endl;
  }
  else {
    MSG("NuFCExperimentFactory",Msg::kInfo) << "Now out of debugging mode.  Systematics are randomly chosen based on input values." << endl;
  }
}


void NuFCExperimentFactory::NoSystematics()
{
  norm_size = 0;
  
  curv_size = 0;
  relcurv_size = 0;
  range_size = 0;
  shower_size = 0;
  relshow_size = 0;
  
  dp_size = 0;
  bg_size = 0;
  
  xsec_size = 0;
  skzp_size = 0;
  
  // Zero out the shifts
  RollSystematics();
}


double NuFCExperimentFactory::GetShift(double size) const
{ 
  if (fixedSyst)     return size;
  else if (gausSyst) return fRandy->Gaus(0, size);
  else               return fRandy->Uniform(-size, size); 
}


void NuFCExperimentFactory::RollSystematics()
{
  // Get a single copy of these systematics
  norm_shift = 1 + GetShift(norm_size);
  
  curv_shift = 1 + GetShift(curv_size);
  relcurv_shift = 1 + GetShift(relcurv_size);
  range_shift = 1 + GetShift(range_size);
  shower_shift = 1 + GetShift(shower_size);
  relshow_shift = 1 + GetShift(relshow_size);

  if (!dp_corr) {
    if (dp_size == 0) dp_shift = 1;
    else if (dp_old)  dp_shift = 1.0 + GetShift(dp_size);
    else              dp_shift = 0.8 + GetShift(dp_size);
    if (dp_shift < 0) dp_shift = 0;
  }
  bg_shift = 1 + GetShift(bg_size);
  
  xsec_shift = GetShift(xsec_size);
  skzp_shift = GetShift(skzp_size);
}


void NuFCExperimentFactory::InvertShifts()
{
  // Get a single copy of these systematics
  norm_shift = 2 - norm_shift;
  
  curv_shift = 2 - curv_shift;
  relcurv_shift = 2 - curv_shift;
  range_shift = 2 - range_shift;
  shower_shift = 2 - shower_shift;
  relshow_shift = 2 - relshow_shift;
  
  dp_shift = 2 - dp_shift;
  bg_shift = 2 - bg_shift;
  
  xsec_shift *= -1;
  skzp_shift *= -1;
}

double NuFCExperimentFactory::GetNormShift() const
{
  return norm_shift;
}

double NuFCExperimentFactory::GetTrkEnergy(const NuFCEvent* nu) const
{
  double E = nu->trkEn;
  if (nu->containedTrk == 1) E *= range_shift;
  else {
    E *= curv_shift;
    if (nu->detector == 2) E *= relcurv_shift;
  }
  return E;
}

double NuFCExperimentFactory::GetShowerEnergy(const NuFCEvent* nu) const
{
  double E = nu->shwEn * shower_shift;
  if (nu->detector == 2) E *= relshow_shift;
  return E;
}


bool NuFCExperimentFactory::IsDP(const NuFCEvent* nu) const
{
  if (nu->ppvz < 4500) return false;
  if (abs(nu->ptype) == 13) return false;
  return true;
  
}

double NuFCExperimentFactory::GetDPShift(const NuFCEvent* nu) const
{
  if (IsDP(nu)) return dp_shift;
  return 1;  
}

double NuFCExperimentFactory::GetBGShift(const NuFCEvent* nu) const
{
  if (nu->iaction == 0 || nu->inu == 14) return bg_shift;
  return 1;
}

double NuFCExperimentFactory::GetBGShift(const Sample_t s) const
{
  if (s == kWrongSignPQ || s == kNCPQ) return bg_shift;
  return 1;
}
      
double NuFCExperimentFactory::GetXSecShift(const NuFCEvent* ev) const
{
  if (ev->iaction != 1) return 1;
  
  static double xp = 25.;
  static double a = 0.895;
  static double xc = 4.0;
  static double b = 7.59e-3;
  static double c = -8.05e-4;
  double E = ev->neuEnMC;
  double shift = -0.0617;
  
  if (E < 25) shift += a - 1 + 2.*(1.-a)*E/xp + (a - 1.)*E*E/(xp*xp);
  if (E < xc) shift += (c*(xc-E)*(xc-E)*(xc-E)+b*(xc-E)*(xc-E));
  
  shift = 1 + xsec_shift * shift;
  return shift;
}

double NuFCExperimentFactory::GetSKZPShift(const NuFCEvent* nu) const
{
  return 1. + skzp_shift*(nu->fluxErr - 1.);
}

---