NuSystematic.cxx

Go to the documentation of this file.

// $Id: NuSystematic.cxx,v 1.28 2008/04/14 15:39:54 evans Exp $
//
// A class to control the systematics studies. Applies systematic
// shifts to NuEvents.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
// $Log: NuSystematic.cxx,v $
// Revision 1.28  2008/04/14 15:39:54  evans
// More unitialised variable fixes for the optimised compiler.
//
// Revision 1.27  2008/03/20 16:52:53  evans
// A few new systematics. Independent ma_res and ma_qe numu scales. A
// combined track range and curvature systematic (assume both are
// correlated, with 2% and 3% shifts respectively; and both are
// furthermore correlated between both detectors). Also an
// energy-dependent form of the shower energy scale systematic.
//
// Revision 1.26  2008/02/19 11:49:57  evans
// Updating the decay pipe parent systematic so it works with the pHE beam.
//
// Revision 1.25  2008/02/17 06:17:00  rhatcher
// Add back some neugen headers that aren't automatically included via
// other headers, was via MCReweight or NeugenWeightCalculator.
//
// Revision 1.24  2008/02/08 16:19:52  evans
// Fixing the BField systematic so it shifts every CC numu event with q/p>0.
//
// Revision 1.23  2008/01/20 01:11:15  evans
// Fixing a bug in the application of the overall antineutrino cross
// section systematic.
//
// Revision 1.22  2008/01/17 01:22:07  evans
// A first attempt at code to systematically shift cross sections.
//
// Revision 1.21  2008/01/14 11:39:16  evans
// Putting in the value I'm going to use for the track fit probability
// systematic.
//
// Revision 1.20  2008/01/14 06:38:30  evans
// A couple of systematics more specific to my selection. One shifts
// track jitter and DPID. The other shifts track fit probability.
//
// Revision 1.19  2008/01/10 18:19:49  evans
// Adding bounds on the Neugen r parameters (resonance/DIS transition
// region) so they can't be shifted below 0 or above 1.
//
// Revision 1.18  2008/01/10 17:32:44  evans
// A number of bug fixes to get the Neugen cross section reweighting
// finally working correctly for every systematic. The shift for the
// numubar resonance-DIS transition region gives a dubiously physical
// reweighting (the -1sigma shift weights some events by -1; the +1sigma
// shift weights them by +3), but it seems to be what's implied by
// DocDB-2989.
//
// Revision 1.17  2008/01/10 01:51:34  evans
// Fixed the Neugen cross section shifting code so it actually returns
// non-unity weights and shifts the Monte Carlo.
//
// Revision 1.16  2008/01/06 16:01:22  evans
// Adding a far-detector-only track energy from curvature shift (4% at
// present; should possibly be 4%+4% in quadrature to be a near-far
// relative shift).
//
// Revision 1.15  2008/01/05 06:23:50  rhatcher
// With new versions of ROOT and gcc (v3.4.3+) we need a #include <cmath>
// in order to have access to sqrt(); perhaps this worked for older versions
// where includes might have been pulled in implicitly.
//
// Revision 1.14  2008/01/04 18:22:48  ahimmel
//
//
// Added in kScraping systematic.  Right now 1 sigma is 25%, but that is just a guess.
//
// I've also added a ShiftAsValue() function analagous to ShiftAsSigma() to reduce a lot of clutter from conversion.  Now the usage of the shift is more obviously indpendent of how it is stored.
//
// Revision 1.13  2008/01/03 16:15:49  evans
// Adding two new systematics: a scale factor in the numubar QEL cross
// sections, and a scale factor in the numubar resonance cross sections
// (both 8%). These replace the numubar CCMA systematic which shouldn't
// be used.
//
// Revision 1.12  2008/01/01 16:46:43  evans
// Adding the ability to shift the NC background by a scale
// factor. Current default is 50%.
//
// Revision 1.11  2007/12/31 18:43:16  evans
// Correcting the separate near and far normalisation code so it applies
// the shift to the correct detector.
//
// Revision 1.10  2007/12/31 12:55:42  evans
// Fixing a bug in the shower energy scale systematics. Adding in
// separate near and far normalisation and B-field systematics.
//
// Revision 1.9  2007/12/30 22:43:09  evans
// Adding a toy B-field systematic (double/half the CC background in the
// nubar selection).
//
// Revision 1.8  2007/12/27 22:11:50  evans
// Debugging and adding a bit of printout.
//
// Revision 1.7  2007/12/23 21:55:41  evans
// Adding in the correct track energy systematics. TrackEnergyScale and
// TrackEnergyOffset no longer exist. Instead, TrackEnergyRange and
// TrackEnergyCurvature are the two systematics to apply (1 sigma = 2%
// and 4% respectively).
//
// Revision 1.6  2007/12/13 22:00:21  evans
// Putting all the cross section systematics in, as detailed in
// DocDB-2989-v6. There are four combined numu+numubar CC cross section
// systematics and three numubar-only CC systematics: so seven systematic
// shifts in all. There are no NC cross section shifts in the class.
//
// Revision 1.5  2007/12/13 14:04:22  evans
// Adding the ability to make a systematic shift of the CCMA
// paramter. All Neugen parameters will follow the same basic method.
//
// Revision 1.4  2007/12/12 16:33:18  evans
// Adding a normalisation systematic of 4% and setting the track energy
// scale systematic to 2%.
//
// Revision 1.3  2007/12/07 15:25:19  evans
// Adding SKZP flux errors. This code takes NuEvent::fluxErr as the
// error, and alters NuEvent::rw with the equation
//
// NuEvent::fluxrw *= 1.0 + sigma*(NuEvent::fluxErr - 1.0).
//
// Revision 1.2  2007/12/06 14:15:36  evans
// Putting in the correct shower energy scale uncertainties.
//
// kShowerEnergyScale is the absolute (near+far simultaneously) scale
// uncertainty. 5.7% from calibration (DocDB-3137), added in quadrature
// with 8% to account for shower modelling uncertainties
// (DocDB-3362). This gives a total 1sigma uncertainty of 10%.
//
// There's also the possibility of doing relative detector
// shifts. kShowerEnergyScaleNear shifts near detector energies by
// 3.1%. kShowerEnergyScaleFar shifts far detector energies by 2.3%. The
// way this systematic is traditionally evaluated is by adding these two
// numbers in quadrature (3.9%) and shifting just one detector (in this
// framework the near detector): kShowerEnergyScaleRelative performs
// this. So either do both the first two of these shifts, or just the
// last one. (All these numbers are calibration-only figures, again from
// DocDB-3137.)
//
// All figures were calculated for the CC box opening in July so may need
// to be updated for the new dataset and for Daikon_04.
//
// Revision 1.1  2007/11/30 17:28:18  evans
// Adding a new class, NuSystematic, which takes a NuEvent and applies
// the relevant sytematic shift to it. At present it only deals with
// shower energy shifts and shower energy offsets (with 1sigma sizes
// picked out of thin air). Takes a NuXMLConfig object in the constructor
// to tell it which systematic and what size/direction of shift to use.
//
//

#include <cmath>
#include "TString.h"

#include "Conventions/BeamType.h"
#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "Conventions/ReleaseType.h"
#include "MessageService/MsgService.h"
#include "MCReweight/MCReweight.h"
#include "MCReweight/NeugenWeightCalculator.h"
#include "NeugenInterface/ccnc.h"
#include "NeugenInterface/flavor.h"
#include "NeugenInterface/init_state.h"
#include "NeugenInterface/interaction.h"
#include "NeugenInterface/nucleus.h"
#include "NeugenInterface/process.h"
#include "NeugenInterface/neugen_cuts.h"
#include "NeugenInterface/neugen_config.h"
#include "NeugenInterface/neugen_wrapper.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuFluctuator.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuXMLConfig.h"
#include "Registry/Registry.h"

ClassImp(NuSystematic)

CVSID("$Id: NuSystematic.cxx,v 1.28 2008/04/14 15:39:54 evans Exp $");

Bool_t NuSystematic::firstMCReweight = true;

//____________________________________________________________________72
NuSystematic::NuSystematic()
{
  MSG("NuSystematic.cxx",Msg::kInfo)
    << "Constructing default NuSystematic object" << endl;

  this->ConfigureSystematics();
  this->ConfigureNeugenDefaults();
  fSystematicID = NuSyst::kUnknown;
  fShiftSizeAsSigma = 0.0;
  fFluctuator = new NuFluctuator();
}

//____________________________________________________________________72
NuSystematic::NuSystematic(const NuXMLConfig& xmlConfig)
{
  this->ConfigureSystematics();
  this->ConfigureNeugenDefaults();
  fFluctuator = 0;

  TString systName = xmlConfig.Name();
  if (systName.Contains("Nominal",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNominal;
  }
  else if (systName.Contains("Scraping",TString::kIgnoreCase) ||
           systName.Contains("DecayPipe",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kScraping;
  }
  else if (systName.Contains("ShowerEnergyOffset",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyOffset;
  }
  else if (systName.Contains("ShowerEnergyScaleBoth",TString::kIgnoreCase) ||
           systName.Contains("AbsoluteEnergyCalibration",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyScale;
  }
  else if (systName.Contains("ShowerEnergyScaleFunctionBoth",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyFunction;
  }
  else if (systName.Contains("ShowerEnergyScaleNear",TString::kIgnoreCase) ||
           systName.Contains("RelativeEnergyCalibrationNear",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyScaleNear;
  }
  else if (systName.Contains("ShowerEnergyScaleFar",TString::kIgnoreCase) ||
           systName.Contains("RelativeEnergyCalibrationFar",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyScaleFar;
  }
  else if (systName.Contains("ShowerEnergyScaleRelative",TString::kIgnoreCase) ||
           systName.Contains("RelativeEnergyCalibration",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kShowerEnergyScaleRelative;
  }
  else if (systName.Contains("TrackEnergyCurvatureBoth",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kTrackEnergyCurvatureBoth;
  }
  else if (systName.Contains("TrackEnergyCurvatureFar",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kTrackEnergyCurvatureFar;
  }
  else if (systName.Contains("TrackEnergyRange",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kTrackEnergyRange;
  }
  else if (systName.Contains("TrackEnergyOverall",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kTrackEnergyOverall;
  }
  else if (systName.Contains("Beam",TString::kIgnoreCase) ||
           systName.Contains("SKZP",TString::kIgnoreCase) ||
           systName.Contains("Flux",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kBeam;
  }
  else if (systName.Contains("NCBackground",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNCBackground;
  }
  else if (systName.Contains("NormalisationBoth",TString::kIgnoreCase) ||
           systName.Contains("NormalizationBoth",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNormalisationBoth;
  }
  else if (systName.Contains("NormalisationNear",TString::kIgnoreCase) ||
           systName.Contains("NormalizationNear",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNormalisationNear;
  }
  else if (systName.Contains("NormalisationFar",TString::kIgnoreCase) ||
           systName.Contains("NormalizationFar",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNormalisationFar;
  }
  else if (systName.Contains("CombinedXSecCCMA",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecCCMA;
  }
  else if (systName.Contains("CombinedXSecMaRes",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecMaRes;
  }
  else if (systName.Contains("CombinedXSecMaQE",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecMaQE;
  }
  else if (systName.Contains("CombinedXSecOverall",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecOverall;
  }
  else if (systName.Contains("CombinedXSecDISMultip2",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecDISMultip2;
  }
  else if (systName.Contains("CombinedXSecDISMultip3",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kCombinedXSecDISMultip3;
  }
  else if (systName.Contains("NuMuBarXSecCCMA",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNuMuBarXSecCCMA;
  }
  else if (systName.Contains("NuMuBarXSecQEL",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNuMuBarXSecQEL;
  }
  else if (systName.Contains("NuMuBarXSecRes",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNuMuBarXSecRes;
  }
  else if (systName.Contains("NuMuBarXSecOverall",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNuMuBarXSecOverall;
  }
  else if (systName.Contains("NuMuBarXSecDISMultip2",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kNuMuBarXSecDISMultip2;
  }
  else if (systName.Contains("BFieldBoth",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kBFieldBoth;
    fFluctuator = new NuFluctuator();
  }
  else if (systName.Contains("BFieldNear",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kBFieldNear;
    fFluctuator = new NuFluctuator();
  }
  else if (systName.Contains("BFieldFar",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kBFieldFar;
    fFluctuator = new NuFluctuator();
  }
  else if (systName.Contains("JitterVDPID",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kJitterVDPID;
  }
  else if (systName.Contains("TFProb",TString::kIgnoreCase)){
    fSystematicID = NuSyst::kTFProb;
  }
  else {
    fSystematicID = NuSyst::kUnknown;
  }

  fShiftSizeAsSigma = this->ConvertShiftToSigma(xmlConfig.Shift());

  MSG("NuSystematic.cxx",Msg::kInfo)
    << "Constructed NuSystematic object for " 
    << xmlConfig.Name() << " "
    << fShiftSizeAsSigma << "sigma" << endl;
}

//____________________________________________________________________72
NuSystematic::~NuSystematic()
{
  if (fFluctuator){delete fFluctuator; fFluctuator = 0;}
}

//____________________________________________________________________72
void NuSystematic::ConfigureSystematics()
{
  fShwEnOff1Sigma = 100*Munits::MeV;
  fShwEnScale1Sigma = 1.10;
  fShwEnFunction1Sigma = 1.0;
  fShwEnScaleFar1Sigma = 1.023;
  fShwEnScaleNear1Sigma = 1.031;
  fShwEnScaleRelative1Sigma = 1.039;
  fTrkEnCurvBoth1Sigma = 1.03;
  fTrkEnCurvFar1Sigma = 1.03;
  fTrkEnRange1Sigma = 1.02;
  fBeam1Sigma = 1.0;
  fBFieldBoth1Sigma = 2.0;
  fBFieldNear1Sigma = 2.0;
  fBFieldFar1Sigma = 2.0;
  fNCBackground1Sigma = 1.50;
  fNormalisationBoth1Sigma = 1.04;
  fNormalisationNear1Sigma = 1.04;
  fNormalisationFar1Sigma = 1.04;
  fCombinedXSecOverall1Sigma = 1.035;
  fCombinedXSecCCMA1Sigma = 1.15;
  fCombinedXSecMaRes1Sigma = 1.15;
  fCombinedXSecMaQE1Sigma = 1.15;
  fCombinedXSecDISMultip21Sigma = 0.1;
  fCombinedXSecDISMultip31Sigma = 0.2;
  fNuMuBarXSecOverall1Sigma = 1.04;
  fNuMuBarXSecCCMA1Sigma = 1.08;
  fNuMuBarXSecQEL1Sigma = 1.08;
  fNuMuBarXSecRes1Sigma = 1.08;
  fNuMuBarXSecDISMultip21Sigma = 0.2;
  fScraping1Sigma = 1.25;
  fJitter1Sigma = 0.01;
  fDPID1Sigma = 0.03;
  fTFProb1Sigma = 0.005;//Complete guess.
}

//____________________________________________________________________72
void NuSystematic::ConfigureNeugenDefaults()
{
  fkno_r112Default = 0.1;
  fkno_r122Default = 0.3;
  fkno_r132Default = 0.3;
  fkno_r142Default = 0.1;
  fkno_r113Default = 1.0;
  fkno_r123Default = 1.0;
  fkno_r133Default = 1.0;
  fkno_r143Default = 1.0;
  fkno_r212Default = 0.1;
  fkno_r222Default = 0.3;
  fkno_r232Default = 0.3;
  fkno_r242Default = 0.1;
  fkno_r213Default = 1.0;
  fkno_r223Default = 1.0;
  fkno_r233Default = 1.0;
  fkno_r243Default = 1.0;
  fma_qeDefault = 0.99;
  fma_resDefault = 1.12;
}

//____________________________________________________________________72
const Float_t NuSystematic::ConvertSigmaToShift(const Float_t sigma)
  const
{
  if (NuSyst::kNominal == fSystematicID){
    return 0.0;
  }
  else if(NuSyst::kScraping == fSystematicID){
    return sigma*(fScraping1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kShowerEnergyOffset == fSystematicID){
    return sigma*fShwEnOff1Sigma;
  }
  else if(NuSyst::kShowerEnergyScale == fSystematicID){
    return sigma*(fShwEnScale1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kShowerEnergyFunction == fSystematicID){
    return sigma;
  }
  else if(NuSyst::kShowerEnergyScaleNear == fSystematicID){
    return sigma*(fShwEnScaleNear1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kShowerEnergyScaleFar == fSystematicID){
    return sigma*(fShwEnScaleFar1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kShowerEnergyScaleRelative == fSystematicID){
    return sigma*(fShwEnScaleRelative1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kTrackEnergyCurvatureBoth == fSystematicID){
    return sigma*(fTrkEnCurvBoth1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kTrackEnergyCurvatureFar == fSystematicID){
    return sigma*(fTrkEnCurvFar1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kTrackEnergyRange == fSystematicID){
    return sigma*(fTrkEnRange1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kTrackEnergyOverall == fSystematicID){
    return sigma;
  }
  else if(NuSyst::kBeam == fSystematicID){
    return sigma;
  }
  else if (NuSyst::kBFieldBoth == fSystematicID){
    return sigma*(fBFieldBoth1Sigma-1.0)+1.0;
  }
  else if (NuSyst::kBFieldNear == fSystematicID){
    return sigma*(fBFieldNear1Sigma-1.0)+1.0;
  }
  else if (NuSyst::kBFieldFar == fSystematicID){
    return sigma*(fBFieldFar1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNCBackground == fSystematicID){
    return sigma*(fNCBackground1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNormalisationBoth == fSystematicID){
    return sigma*(fNormalisationBoth1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNormalisationNear == fSystematicID){
    return sigma*(fNormalisationNear1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNormalisationFar == fSystematicID){
    return sigma*(fNormalisationFar1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kCombinedXSecOverall == fSystematicID){
    return sigma*(fCombinedXSecOverall1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kCombinedXSecCCMA == fSystematicID){
    return sigma*(fCombinedXSecCCMA1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kCombinedXSecMaRes == fSystematicID){
    return sigma*(fCombinedXSecMaRes1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kCombinedXSecMaQE == fSystematicID){
    return sigma*(fCombinedXSecMaQE1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kCombinedXSecDISMultip2 == fSystematicID){
    return sigma*fCombinedXSecDISMultip21Sigma;
  }
  else if(NuSyst::kCombinedXSecDISMultip3 == fSystematicID){
    return sigma*fCombinedXSecDISMultip31Sigma;
  }
  else if(NuSyst::kNuMuBarXSecOverall == fSystematicID){
    return sigma*(fNuMuBarXSecOverall1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNuMuBarXSecCCMA == fSystematicID){
    return sigma*(fNuMuBarXSecCCMA1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNuMuBarXSecQEL == fSystematicID){
    return sigma*(fNuMuBarXSecQEL1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNuMuBarXSecRes == fSystematicID){
    return sigma*(fNuMuBarXSecRes1Sigma-1.0)+1.0;
  }
  else if(NuSyst::kNuMuBarXSecDISMultip2 == fSystematicID){
    return sigma*fNuMuBarXSecDISMultip21Sigma;
  }
  else if(NuSyst::kJitterVDPID == fSystematicID){
    return sigma;
  }
  else if(NuSyst::kTFProb == fSystematicID){
    return sigma;
  }
  else {
    return 0.0;
  }
}

//____________________________________________________________________72
const Float_t NuSystematic::ConvertShiftToSigma(const Float_t shift)
  const
{
  if (NuSyst::kNominal == fSystematicID){
    return 0.0;
  }
  else if(NuSyst::kScraping == fSystematicID){
    return (shift-1.0)/(fScraping1Sigma-1.0);
  }
  else if(NuSyst::kShowerEnergyOffset == fSystematicID){
    return shift/fShwEnOff1Sigma;
  }
  else if(NuSyst::kShowerEnergyScale == fSystematicID){
    return (shift-1.0)/(fShwEnScale1Sigma-1.0);
  }
  else if(NuSyst::kShowerEnergyFunction == fSystematicID){
    return shift;
  }
  else if(NuSyst::kShowerEnergyScaleNear == fSystematicID){
    return (shift-1.0)/(fShwEnScaleNear1Sigma-1.0);
  }
  else if(NuSyst::kShowerEnergyScaleFar == fSystematicID){
    return (shift-1.0)/(fShwEnScaleFar1Sigma-1.0);
  }
  else if(NuSyst::kShowerEnergyScaleRelative == fSystematicID){
    return (shift-1.0)/(fShwEnScaleRelative1Sigma-1.0);
  }
  else if(NuSyst::kTrackEnergyRange == fSystematicID){
    return (shift-1.0)/(fTrkEnRange1Sigma-1.0);
  }
  else if(NuSyst::kTrackEnergyCurvatureBoth == fSystematicID){
    return (shift-1.0)/(fTrkEnCurvBoth1Sigma-1.0);
  }
  else if(NuSyst::kTrackEnergyCurvatureFar == fSystematicID){
    return (shift-1.0)/(fTrkEnCurvFar1Sigma-1.0);
  }
  else if(NuSyst::kTrackEnergyOverall == fSystematicID){
    return shift;
  }
  else if(NuSyst::kBeam == fSystematicID){
    return shift;
  }
  else if (NuSyst::kBFieldBoth == fSystematicID){
    return (shift-1.0)/(fBFieldBoth1Sigma-1.0);
  }
  else if (NuSyst::kBFieldNear == fSystematicID){
    return (shift-1.0)/(fBFieldNear1Sigma-1.0);
  }
  else if (NuSyst::kBFieldFar == fSystematicID){
    return (shift-1.0)/(fBFieldFar1Sigma-1.0);
  }
  else if(NuSyst::kNCBackground == fSystematicID){
    return (shift-1.0)/(fNCBackground1Sigma-1.0);
  }
  else if(NuSyst::kNormalisationBoth == fSystematicID){
    return (shift-1.0)/(fNormalisationBoth1Sigma-1.0);
  }
  else if(NuSyst::kNormalisationNear == fSystematicID){
    return (shift-1.0)/(fNormalisationNear1Sigma-1.0);
  }
  else if(NuSyst::kNormalisationFar == fSystematicID){
    return (shift-1.0)/(fNormalisationFar1Sigma-1.0);
  }
  else if(NuSyst::kCombinedXSecOverall == fSystematicID){
    return (shift-1.0)/(fCombinedXSecOverall1Sigma-1.0);
  }
  else if(NuSyst::kCombinedXSecCCMA == fSystematicID){
    return (shift-1.0)/(fCombinedXSecCCMA1Sigma-1.0);
  }
  else if(NuSyst::kCombinedXSecMaRes == fSystematicID){
    return (shift-1.0)/(fCombinedXSecMaRes1Sigma-1.0);
  }
  else if(NuSyst::kCombinedXSecMaQE == fSystematicID){
    return (shift-1.0)/(fCombinedXSecMaQE1Sigma-1.0);
  }
  else if(NuSyst::kCombinedXSecDISMultip2 == fSystematicID){
    return shift/fCombinedXSecDISMultip21Sigma;
  }
  else if(NuSyst::kCombinedXSecDISMultip3 == fSystematicID){
    return shift/fCombinedXSecDISMultip31Sigma;
  }
  else if(NuSyst::kNuMuBarXSecOverall == fSystematicID){
    return (shift-1.0)/(fNuMuBarXSecOverall1Sigma-1.0);
  }
  else if(NuSyst::kNuMuBarXSecQEL == fSystematicID){
    return (shift-1.0)/(fNuMuBarXSecQEL1Sigma-1.0);
  }
  else if(NuSyst::kNuMuBarXSecRes == fSystematicID){
    return (shift-1.0)/(fNuMuBarXSecRes1Sigma-1.0);
  }
  else if(NuSyst::kNuMuBarXSecDISMultip2 == fSystematicID){
    return shift/fNuMuBarXSecDISMultip21Sigma;
  }
  else if(NuSyst::kJitterVDPID == fSystematicID){
    return shift;
  }
  else if(NuSyst::kTFProb == fSystematicID){
    return shift;
  }
  else {
    return 0.0;
  }
}

//____________________________________________________________________72
void NuSystematic::Shift(NuEvent& event) const
{
  if (NuSyst::kNominal == fSystematicID){
    return;
  }
  else if(NuSyst::kScraping == fSystematicID){
    this->ScrapingShift(event);
    return;
  }
  else if(NuSyst::kShowerEnergyOffset == fSystematicID){
    this->ShowerEnergyOffset(event);
    return;
  }
  else if(NuSyst::kShowerEnergyScale == fSystematicID){
    this->ShowerEnergyScale(event);
    return;
  }
  else if(NuSyst::kShowerEnergyFunction == fSystematicID){
    this->ShowerEnergyFunction(event);
    return;
  }
  else if(NuSyst::kShowerEnergyScaleNear == fSystematicID){
    if (Detector::kNear==event.detector){
      this->ShowerEnergyScale(event);
    }
    return;
  }
  else if(NuSyst::kShowerEnergyScaleFar == fSystematicID){
    if (Detector::kFar==event.detector){
      this->ShowerEnergyScale(event);
    }
    return;
  }
  else if(NuSyst::kShowerEnergyScaleRelative == fSystematicID){
    if (Detector::kNear==event.detector){
      this->ShowerEnergyScale(event);
    }
    return;
  }
  else if(NuSyst::kTrackEnergyRange == fSystematicID){
    if (event.usedRange){
      this->TrackEnergyRange(event);
    }
    return;
  }
  else if(NuSyst::kTrackEnergyCurvatureBoth == fSystematicID){
    if (event.usedCurv){
      this->TrackEnergyCurvatureBoth(event);
    }
    return;
  }
  else if(NuSyst::kTrackEnergyCurvatureFar == fSystematicID){
    if (event.usedCurv && Detector::kFar==event.detector){
      this->TrackEnergyCurvatureFar(event);
    }
    return;
  }
  else if(NuSyst::kTrackEnergyOverall == fSystematicID){
    this->TrackEnergyOverall(event);
    return;
  }
  else if(NuSyst::kBeam == fSystematicID){
    this->BeamShift(event);
    return;
  }
  else if (NuSyst::kBFieldBoth == fSystematicID){
    this->BFieldBothShift(event);
  }
  else if (NuSyst::kBFieldNear == fSystematicID){
    if (Detector::kNear==event.detector){
      this->BFieldNearShift(event);
    }
    return;
  }
  else if (NuSyst::kBFieldFar == fSystematicID){
    if (Detector::kFar==event.detector){
      this->BFieldFarShift(event);
    }
    return;
  }
  else if(NuSyst::kNCBackground == fSystematicID){
    this->NCBackgroundShift(event);
    return;
  }
  else if(NuSyst::kNormalisationBoth == fSystematicID){
    this->NormalisationBothShift(event);
    return;
  }
  else if(NuSyst::kNormalisationNear == fSystematicID){
    if (Detector::kNear == event.detector){
      this->NormalisationNearShift(event);
    }
    return;
  }
  else if(NuSyst::kNormalisationFar == fSystematicID){
    if (Detector::kFar == event.detector){
      this->NormalisationFarShift(event);
    }
    return;
  }
  else if(NuSyst::kCombinedXSecOverall == fSystematicID ||
          NuSyst::kNuMuBarXSecOverall == fSystematicID){
    this->OverallXSecShift(event);
    return;
  }
  else if(NuSyst::kNuMuBarXSecQEL == fSystematicID ||
          NuSyst::kNuMuBarXSecQEL == fSystematicID){
    this->NuMuBarQELXSecShift(event);
    return;
  }
  else if(NuSyst::kNuMuBarXSecRes == fSystematicID ||
          NuSyst::kNuMuBarXSecRes == fSystematicID){
    this->NuMuBarResXSecShift(event);
    return;
  }
  else if(NuSyst::kCombinedXSecCCMA == fSystematicID ||
          NuSyst::kCombinedXSecMaRes == fSystematicID ||
          NuSyst::kCombinedXSecMaQE == fSystematicID ||
          NuSyst::kCombinedXSecDISMultip2 == fSystematicID ||
          NuSyst::kCombinedXSecDISMultip3 == fSystematicID ||
          NuSyst::kNuMuBarXSecCCMA == fSystematicID ||
          NuSyst::kNuMuBarXSecDISMultip2 == fSystematicID){
    this->NeugenXSecShift(event);
    return;
  }
  else if(NuSyst::kJitterVDPID == fSystematicID){
    this->JitterVDPIDShift(event);
    return;
  }
  else {
    return;
  }
}

//____________________________________________________________________72
void NuSystematic::BeamShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing beam shift" << endl;
  event.rw *= 1.0 +
    this->ConvertSigmaToShift(this->ShiftAsSigma()) *
    (event.fluxErr-1.0);
  return;
}

//____________________________________________________________________72
void NuSystematic::JitterVDPIDShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track jitter & DPID shift" << endl;
  event.jitter += this->ShiftAsSigma()*fJitter1Sigma;
  event.dpID += this->ShiftAsSigma()*fDPID1Sigma;
  return;
}

//____________________________________________________________________72
void NuSystematic::TFProbShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track fit probability shift" << endl;
  event.prob += this->ShiftAsSigma()*fTFProb1Sigma;
  return;
}

//____________________________________________________________________72
void NuSystematic::NeugenXSecShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing a Neugen cross section shift" << endl;
  if (0 == event.iaction){return;}
  if (NuSyst::kNuMuBarXSecCCMA == this->SystematicID() ||
      NuSyst::kNuMuBarXSecDISMultip2 == this->SystematicID()){
    if (14 == event.inu){return;}
  }

  MCReweight& mcReweight = MCReweight::Instance();
  if (firstMCReweight){
    NeugenWeightCalculator* wc = new NeugenWeightCalculator();
    mcReweight.AddWeightCalculator(wc);
    firstMCReweight = false;
  }
  //  cout << "Num weight calc: " << mcReweight.NumWeightCalcAdded() << endl;
  Registry reweightConfigRegistry;
  Registry defaultRegistry;

  //Set to MODBYRS4 for Daikon.
  if (ReleaseType::IsDaikon(event.releaseType)){
    reweightConfigRegistry.Set("neugen:config_name","MODBYRS");
    reweightConfigRegistry.Set("neugen:config_no",4);
    defaultRegistry.Set("neugen:config_name","MODBYRS");
    defaultRegistry.Set("neugen:config_no",4);
  }
  else{
    MSG("NuSystematic.cxx",Msg::kError)
      << "Using non-daikon MC. I don't know how to apply Neugen "
      << "parameters to that."
      << endl;
  }

  //Set the shifted CombinedXSecCCMA parameters.
  this->SetShiftedNeugenParameters(reweightConfigRegistry,
                                   this->SystematicID());

  //Create objects to give to MCReweight.
  MCEventInfo mcEventInfo = this->CreateMCEventInfo(event);
  NuParent* nuParent = 0;

  mcReweight.ResetAllReweightConfigs();

  //    cout << "iresonance is " << event.iresonance << endl;
  //Get the weight.
  Double_t weight = mcReweight.ComputeWeight(&mcEventInfo,
                                             nuParent,
                                             &reweightConfigRegistry);
  //  mcReweight.PrintReweightConfig(cout);


  //Get the default weight.
  this->SetNeugenDefaults(defaultRegistry);
  mcReweight.ResetAllReweightConfigs();
  Double_t defaultWeight = mcReweight.ComputeWeight(&mcEventInfo,
                                                    nuParent,
                                                    &defaultRegistry);
  //  cout << "Weight changed by " << weight << "/" << defaultWeight << endl;
  if (defaultWeight>0.0){
    event.rw *= weight/defaultWeight;
  }
  else{
    MSG("NuSystematic.cxx",Msg::kError)
      << "Default weight <= 0."
      << endl;
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::OverallXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing overall cross section shift" << endl;
  if (1 != event.iaction){return;}
  if (!((14 == event.inu) || (-14 == event.inu))){return;}
  if (NuSyst::kNuMuBarXSecOverall == this->SystematicID()){
    if (-14 != event.inu){return;}
  }
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::NuMuBarQELXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NuMuBar QEL cross section shift" << endl;
  if (1 != event.iaction){return;}
  if (-14 != event.inu){return;}
  if (1001 != event.iresonance){return;}
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::NuMuBarResXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NuMuBar Res cross section shift" << endl;
  if (1 != event.iaction){return;}
  if (-14 != event.inu){return;}
  if (1002 != event.iresonance){return;}
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::BFieldBothShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing BFieldBoth shift"
    << this->ShiftAsValue() << endl;
  if (1==event.charge
      && 14 == event.inu && 1 == event.iaction){
    event.rw *= this->ShiftAsValue();
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::BFieldNearShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing BFieldNear shift"
    << this->ShiftAsValue() << endl;
  if (1==event.charge
      && 14 == event.inu && 1 == event.iaction){
    event.rw *= this->ShiftAsValue();
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::BFieldFarShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,10)
    << "Performing BFieldFar shift "
    << this->ShiftAsValue() << endl;
//   if (NuFluctuator::kCCNuMuBar == fFluctuator->SelectedAs(event)
  if (1==event.charge
      && 14 == event.inu && 1 == event.iaction){
    event.rw *= this->ShiftAsValue();
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::NCBackgroundShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NC background shift" << endl;
  if (0 != event.iaction){return;}
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::NormalisationBothShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing NormalisationBoth shift" << endl;
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::NormalisationNearShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing NormalisationNear shift" << endl;
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::NormalisationFarShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing NormalisationFar shift" << endl;
  event.rw *= this->ShiftAsValue();
  return;
}

//____________________________________________________________________72
void NuSystematic::ScrapingShift(NuEvent& event) const 
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing scraping shift" << endl;
    Float_t r = sqrt(event.ppvx*event.ppvx +
                     event.ppvy*event.ppvy);
    Float_t z = event.ppvz;
    Float_t Znom = 52.06; // -187.06 for HE
    if (BeamType::kL250z200i == event.beamType){
      Znom = -187.06;
    }
    else if (BeamType::kL010z185i == event.beamType){
      Znom = 52.06;
    }
    else{
      Znom = 52.06;
    }
    
    if (r < 1.65 && TMath::Abs(z - Znom) < 0.1)
        return; // Target End
    
    if (TMath::Abs(r - 1.51) < 0.11 && z < Znom)
        return; // Target Side
    
    // Chase (z < 4500) or Decay Pipe (z > 4500)
    event.rw *= this->ShiftAsValue();
    
    return;
}

//____________________________________________________________________72
void NuSystematic::ShowerEnergyOffset(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing shower energy offset shift" << endl;
  event.shwEn +=
    this->ConvertSigmaToShift(this->ShiftAsSigma())/Munits::GeV;
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::ShowerEnergyFunction(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing shower energy function scale " 
    << this->ConvertSigmaToShift(this->ShiftAsSigma())
    << " sigma" << endl;

  Double_t offset = 7.0;
  Double_t scale = 4.0;
  Double_t eLife = 1.5; //GeV
  Double_t enShift = 1.0 + this->ShiftAsSigma()*0.01*
    (offset + scale*TMath::Exp(-1.0*event.shwEn/eLife));
  event.shwEn *= enShift;
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::ShowerEnergyScale(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing shower energy scale shift" << endl;
  event.shwEn *= this->ShiftAsValue();
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::TrackEnergyRange(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing track energy range shift" << endl;
  event.trkEn *= this->ShiftAsValue();
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::TrackEnergyCurvatureBoth(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing track energy curvature both shift" << endl;
  event.trkEn *= this->ShiftAsValue();
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::TrackEnergyCurvatureFar(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing track energy curvature far shift" << endl;
  event.trkEn *= this->ShiftAsValue();
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::TrackEnergyOverall(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing overall track energy shift "
    << this->ShiftAsSigma() << " sigma"
    << endl;
  Double_t eScale = 1.0;
  if (event.usedRange){
    eScale = 2.0;
  }
  else if (event.usedCurv){
    eScale = 3.0;
  }
  else {
    eScale = 1.0;
  }
  event.trkEn *= 1.0 + this->ShiftAsSigma()*0.01*eScale;
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::SetNeugenDefaults(Registry& registry) const
{
  registry.Set("neugen:kno_r112",fkno_r112Default);
  registry.Set("neugen:kno_r122",fkno_r122Default);
  registry.Set("neugen:kno_r132",fkno_r132Default);
  registry.Set("neugen:kno_r142",fkno_r142Default);
  registry.Set("neugen:kno_r113",fkno_r113Default);
  registry.Set("neugen:kno_r123",fkno_r123Default);
  registry.Set("neugen:kno_r133",fkno_r133Default);
  registry.Set("neugen:kno_r143",fkno_r143Default);
  registry.Set("neugen:ma_qe",fma_qeDefault);
  registry.Set("neugen:ma_res",fma_resDefault);
  return;
}

//____________________________________________________________________72
void NuSystematic
::SetShiftedNeugenParameters(Registry& registry,
                             const NuSyst::NuSystematic_t syst) const
{
  if (NuSyst::kCombinedXSecCCMA == syst){
//     registry.Set("neugen:kno_r112",this->KNO112Default());
//     registry.Set("neugen:kno_r122",this->KNO122Default());
//     registry.Set("neugen:kno_r132",this->KNO132Default());
//     registry.Set("neugen:kno_r142",this->KNO142Default());
//     registry.Set("neugen:kno_r113",this->KNO113Default());
//     registry.Set("neugen:kno_r123",this->KNO123Default());
//     registry.Set("neugen:kno_r133",this->KNO133Default());
//     registry.Set("neugen:kno_r143",this->KNO143Default());
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    registry.Set("neugen:ma_qe",ma_qe);
    registry.Set("neugen:ma_res",ma_res);
    return;
  }
  else if (NuSyst::kCombinedXSecMaRes == syst){
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    registry.Set("neugen:ma_res",ma_res);
    return;
  }
  else if (NuSyst::kCombinedXSecMaQE == syst){
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    registry.Set("neugen:ma_qe",ma_qe);
    return;
  }
  else if (NuSyst::kCombinedXSecDISMultip2 == syst){
//     registry.Set("neugen:ma_qe",this->MA_QEDefault());
//     registry.Set("neugen:ma_res",this->MA_ResDefault());
//     registry.Set("neugen:kno_r113",this->KNO113Default());
//     registry.Set("neugen:kno_r123",this->KNO123Default());
//     registry.Set("neugen:kno_r133",this->KNO133Default());
//     registry.Set("neugen:kno_r143",this->KNO143Default());
    Float_t kno112 = this->KNO112Default() + this->ShiftAsValue();
    if (kno112<0.0){kno112=0.0;}
    if (kno112>1.0){kno112=1.0;}
    Float_t kno122 = this->KNO122Default() + this->ShiftAsValue();
    if (kno122<0.0){kno122=0.0;}
    if (kno122>1.0){kno122=1.0;}
    Float_t kno132 = this->KNO132Default() + this->ShiftAsValue();
    if (kno132<0.0){kno132=0.0;}
    if (kno132>1.0){kno132=1.0;}
    Float_t kno142 = this->KNO142Default() + this->ShiftAsValue();
    if (kno142<0.0){kno142=0.0;}
    if (kno142>1.0){kno142=1.0;}
    Float_t kno212 = this->KNO212Default() + this->ShiftAsValue();
    if (kno212<0.0){kno212=0.0;}
    if (kno212>1.0){kno212=1.0;}
    Float_t kno222 = this->KNO222Default() + this->ShiftAsValue();
    if (kno222<0.0){kno222=0.0;}
    if (kno222>1.0){kno222=1.0;}
    Float_t kno232 = this->KNO232Default() + this->ShiftAsValue();
    if (kno232<0.0){kno232=0.0;}
    if (kno232>1.0){kno232=1.0;}
    Float_t kno242 = this->KNO242Default() + this->ShiftAsValue();
    if (kno242<0.0){kno242=0.0;}
    if (kno242>1.0){kno242=1.0;}
    registry.Set("neugen:kno_r112",kno112);
    registry.Set("neugen:kno_r122",kno122);
    registry.Set("neugen:kno_r132",kno132);
    registry.Set("neugen:kno_r142",kno142);
    registry.Set("neugen:kno_r212",kno212);
    registry.Set("neugen:kno_r222",kno222);
    registry.Set("neugen:kno_r232",kno232);
    registry.Set("neugen:kno_r242",kno242);
  }
  else if (NuSyst::kCombinedXSecDISMultip3 == syst){
//     registry.Set("neugen:ma_qe",this->MA_QEDefault());
//     registry.Set("neugen:ma_res",this->MA_ResDefault());
//     registry.Set("neugen:kno_r112",this->KNO112Default());
//     registry.Set("neugen:kno_r122",this->KNO122Default());
//     registry.Set("neugen:kno_r132",this->KNO132Default());
//     registry.Set("neugen:kno_r142",this->KNO142Default());
    Float_t kno113 = this->KNO113Default() + this->ShiftAsValue();
    if (kno113<0.0){kno113=0.0;}
    if (kno113>1.0){kno113=1.0;}
    Float_t kno123 = this->KNO123Default() + this->ShiftAsValue();
    if (kno123<0.0){kno123=0.0;}
    if (kno123>1.0){kno123=1.0;}
    Float_t kno133 = this->KNO133Default() + this->ShiftAsValue();
    if (kno133<0.0){kno133=0.0;}
    if (kno133>1.0){kno133=1.0;}
    Float_t kno143 = this->KNO143Default() + this->ShiftAsValue();
    if (kno143<0.0){kno143=0.0;}
    if (kno143>1.0){kno143=1.0;}
    Float_t kno213 = this->KNO213Default() + this->ShiftAsValue();
    if (kno213<0.0){kno213=0.0;}
    if (kno213>1.0){kno213=1.0;}
    Float_t kno223 = this->KNO223Default() + this->ShiftAsValue();
    if (kno223<0.0){kno223=0.0;}
    if (kno223>1.0){kno223=1.0;}
    Float_t kno233 = this->KNO233Default() + this->ShiftAsValue();
    if (kno233<0.0){kno233=0.0;}
    if (kno233>1.0){kno233=1.0;}
    Float_t kno243 = this->KNO243Default() + this->ShiftAsValue();
    if (kno243<0.0){kno243=0.0;}
    if (kno243>1.0){kno243=1.0;}
    registry.Set("neugen:kno_r113",kno113);
    registry.Set("neugen:kno_r123",kno123);
    registry.Set("neugen:kno_r133",kno133);
    registry.Set("neugen:kno_r143",kno143);
    registry.Set("neugen:kno_r213",kno213);
    registry.Set("neugen:kno_r223",kno223);
    registry.Set("neugen:kno_r233",kno233);
    registry.Set("neugen:kno_r243",kno243);
  }
  else if (NuSyst::kNuMuBarXSecCCMA == syst){
//     registry.Set("neugen:kno_r112",this->KNO112Default());
//     registry.Set("neugen:kno_r122",this->KNO122Default());
//     registry.Set("neugen:kno_r132",this->KNO132Default());
//     registry.Set("neugen:kno_r142",this->KNO142Default());
//     registry.Set("neugen:kno_r113",this->KNO113Default());
//     registry.Set("neugen:kno_r123",this->KNO123Default());
//     registry.Set("neugen:kno_r133",this->KNO133Default());
//     registry.Set("neugen:kno_r143",this->KNO143Default());
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    registry.Set("neugen:ma_qe",ma_qe);
    registry.Set("neugen:ma_res",ma_res);
    return;
  }
  else if (NuSyst::kNuMuBarXSecDISMultip2 == syst){
//     registry.Set("neugen:ma_qe",this->MA_QEDefault());
//     registry.Set("neugen:ma_res",this->MA_ResDefault());
//     registry.Set("neugen:kno_r112",this->KNO112Default());
//     registry.Set("neugen:kno_r122",this->KNO122Default());
//     registry.Set("neugen:kno_r113",this->KNO113Default());
//     registry.Set("neugen:kno_r123",this->KNO123Default());
//     registry.Set("neugen:kno_r133",this->KNO133Default());
//     registry.Set("neugen:kno_r143",this->KNO143Default());
    Float_t kno132 = this->KNO132Default() + this->ShiftAsValue();
    if (kno132<0.0){kno132=0.0;}
    if (kno132>1.0){kno132=1.0;}
    Float_t kno142 = this->KNO142Default() + this->ShiftAsValue();
    if (kno142<0.0){kno142=0.0;}
    if (kno142>1.0){kno142=1.0;}
    registry.Set("neugen:kno_r132",kno132);
    registry.Set("neugen:kno_r142",kno142);
  }
}

//____________________________________________________________________72
const MCEventInfo NuSystematic
::CreateMCEventInfo(const NuEvent& nuEvent) const
{
  MCEventInfo mcEvent;

  mcEvent.nuE = nuEvent.neuEnMC;
  mcEvent.nuPx = nuEvent.neuPxMC;
  mcEvent.nuPy = nuEvent.neuPyMC;
  mcEvent.nuPz = nuEvent.neuPzMC;
  mcEvent.tarE = nuEvent.tgtEnMC;
  mcEvent.tarPx = nuEvent.tgtPxMC;
  mcEvent.tarPy = nuEvent.tgtPyMC;
  mcEvent.tarPz = nuEvent.tgtPzMC;
  mcEvent.y = nuEvent.yMC;
  mcEvent.x = nuEvent.xMC;
  mcEvent.q2 = nuEvent.q2MC;
  mcEvent.w2 = nuEvent.w2MC;
  mcEvent.iaction = nuEvent.iaction;
  mcEvent.inu = nuEvent.inu;
  mcEvent.iresonance = nuEvent.iresonance;
  mcEvent.initial_state = nuEvent.initialStateMC;
  mcEvent.nucleus = nuEvent.nucleusMC;
  mcEvent.had_fs = nuEvent.hadronicFinalStateMC;

  return mcEvent;
}

//____________________________________________________________________72
void NuSystematic
::SetShiftedNeugenParameters(neugen_config& config) const
{
  NuSyst::NuSystematic_t syst = this->SystematicID();
  if (NuSyst::kCombinedXSecCCMA == syst){
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    config.set_ma_qe(ma_qe);
    config.set_ma_res(ma_res);
    return;
  }
  else if (NuSyst::kCombinedXSecMaRes == syst){
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    config.set_ma_res(ma_res);
    return;
  }
  else if (NuSyst::kCombinedXSecMaQE == syst){
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    config.set_ma_qe(ma_qe);
    return;
  }
  else if (NuSyst::kCombinedXSecDISMultip2 == syst){
    Float_t kno112 = this->KNO112Default() + this->ShiftAsValue();
    if (kno112<0.0){kno112=0.0;}
    if (kno112>1.0){kno112=1.0;}
    Float_t kno122 = this->KNO122Default() + this->ShiftAsValue();
    if (kno122<0.0){kno122=0.0;}
    if (kno122>1.0){kno122=1.0;}
    Float_t kno132 = this->KNO132Default() + this->ShiftAsValue();
    if (kno132<0.0){kno132=0.0;}
    if (kno132>1.0){kno132=1.0;}
    Float_t kno142 = this->KNO142Default() + this->ShiftAsValue();
    if (kno142<0.0){kno142=0.0;}
    if (kno142>1.0){kno142=1.0;}
    Float_t kno212 = this->KNO212Default() + this->ShiftAsValue();
    if (kno212<0.0){kno212=0.0;}
    if (kno212>1.0){kno212=1.0;}
    Float_t kno222 = this->KNO222Default() + this->ShiftAsValue();
    if (kno222<0.0){kno222=0.0;}
    if (kno222>1.0){kno222=1.0;}
    Float_t kno232 = this->KNO232Default() + this->ShiftAsValue();
    if (kno232<0.0){kno232=0.0;}
    if (kno232>1.0){kno232=1.0;}
    Float_t kno242 = this->KNO242Default() + this->ShiftAsValue();
    if (kno242<0.0){kno242=0.0;}
    if (kno242>1.0){kno242=1.0;}
    config.set_dis_res(1,2,e_vp,kno112);
    config.set_dis_res(1,2,e_vn,kno122);
    config.set_dis_res(1,2,e_vbp,kno132);
    config.set_dis_res(1,2,e_vbn,kno142);
    config.set_dis_res(2,2,e_vp,kno212);
    config.set_dis_res(2,2,e_vn,kno222);
    config.set_dis_res(2,2,e_vbp,kno232);
    config.set_dis_res(2,2,e_vbn,kno242);
  }
  else if (NuSyst::kCombinedXSecDISMultip3 == syst){
    Float_t kno113 = this->KNO113Default() + this->ShiftAsValue();
    if (kno113<0.0){kno113=0.0;}
    if (kno113>1.0){kno113=1.0;}
    Float_t kno123 = this->KNO123Default() + this->ShiftAsValue();
    if (kno123<0.0){kno123=0.0;}
    if (kno123>1.0){kno123=1.0;}
    Float_t kno133 = this->KNO133Default() + this->ShiftAsValue();
    if (kno133<0.0){kno133=0.0;}
    if (kno133>1.0){kno133=1.0;}
    Float_t kno143 = this->KNO143Default() + this->ShiftAsValue();
    if (kno143<0.0){kno143=0.0;}
    if (kno143>1.0){kno143=1.0;}
    Float_t kno213 = this->KNO213Default() + this->ShiftAsValue();
    if (kno213<0.0){kno213=0.0;}
    if (kno213>1.0){kno213=1.0;}
    Float_t kno223 = this->KNO223Default() + this->ShiftAsValue();
    if (kno223<0.0){kno223=0.0;}
    if (kno223>1.0){kno223=1.0;}
    Float_t kno233 = this->KNO233Default() + this->ShiftAsValue();
    if (kno233<0.0){kno233=0.0;}
    if (kno233>1.0){kno233=1.0;}
    Float_t kno243 = this->KNO243Default() + this->ShiftAsValue();
    if (kno243<0.0){kno243=0.0;}
    if (kno243>1.0){kno243=1.0;}
    config.set_dis_res(1,3,e_vp,kno113);
    config.set_dis_res(1,3,e_vn,kno123);
    config.set_dis_res(1,3,e_vbp,kno133);
    config.set_dis_res(1,3,e_vbn,kno143);
    config.set_dis_res(2,3,e_vp,kno213);
    config.set_dis_res(2,3,e_vn,kno223);
    config.set_dis_res(2,3,e_vbp,kno233);
    config.set_dis_res(2,3,e_vbn,kno243);
  }
  else if (NuSyst::kNuMuBarXSecCCMA == syst){
    Float_t ma_qe = this->MA_QEDefault() * this->ShiftAsValue();
    Float_t ma_res = this->MA_ResDefault() * this->ShiftAsValue();
    config.set_ma_qe(ma_qe);
    config.set_ma_res(ma_res);
    return;
  }
  else if (NuSyst::kNuMuBarXSecDISMultip2 == syst){
    Float_t kno132 = this->KNO132Default() + this->ShiftAsValue();
    if (kno132<0.0){kno132=0.0;}
    if (kno132>1.0){kno132=1.0;}
    Float_t kno142 = this->KNO142Default() + this->ShiftAsValue();
    if (kno142<0.0){kno142=0.0;}
    if (kno142>1.0){kno142=1.0;}
    config.set_dis_res(1,2,e_vbp,kno132);
    config.set_dis_res(1,2,e_vbn,kno142);
  }
}

//____________________________________________________________________72
const Double_t NuSystematic
::XSecShiftScale(const Double_t energy,
                 const NuParticle::NuParticleType_t particle) const
{
  if (!(NuParticle::kNuMu == particle || NuParticle::kNuMuBar == particle)){
    MSG("NuSystematic.cxx",Msg::kWarning)
      << "I don't know how to systematically shift this particle"
      << endl;
    return 1.0;
  }

  if (NuSyst::kCombinedXSecOverall == this->SystematicID()){
    return this->ShiftAsValue();
  }
  else if (NuSyst::kNuMuBarXSecOverall == this->SystematicID()){
    if (NuParticle::kNuMuBar == particle){
      return this->ShiftAsValue();
    }
    else {
      return 1.0;
    }
  }
  else if (NuSyst::kNuMuBarXSecCCMA == this->SystematicID() ||
           NuSyst::kNuMuBarXSecDISMultip2 == this->SystematicID()){
    if (NuParticle::kNuMuBar != particle){
      return 1.0;
    }
    else {return this->NeugenXSecShiftScale(energy,particle);}
  }
  else if(NuSyst::kCombinedXSecCCMA == this->SystematicID() ||
          NuSyst::kCombinedXSecDISMultip2 == this->SystematicID() ||
          NuSyst::kCombinedXSecDISMultip3 == this->SystematicID()){
    if (!(NuParticle::kNuMu == particle || NuParticle::kNuMuBar == particle)){
      return 1.0;
    }
    else if(NuSyst::kNuMuBarXSecQEL == this->SystematicID()){
      return this->QELXSecShiftScale(energy,particle);
    }
    else if (NuSyst::kNuMuBarXSecRes == this->SystematicID()){
      return this->ResXSecShiftScale(energy,particle);
    }
    else {return this->NeugenXSecShiftScale(energy,particle);}
  }
  else {return 1.0;}
}

//____________________________________________________________________72
const Double_t NuSystematic
::QELXSecShiftScale(const Double_t energy,
                    const NuParticle::NuParticleType_t particle)
  const 
{
  neugen_config defaultConfig;
  neugen_wrapper defaultWrapper(&defaultConfig);
  init_state_t inlStateP;
  init_state_t inlStateN;
  if (NuParticle::kNuMu == particle){
    inlStateP = e_vp;
    inlStateN = e_vn;
  }
  else if (NuParticle::kNuMuBar == particle){
    inlStateP = e_vbp;
    inlStateN = e_vbn;
  }
  else{
    MSG("NuSystematic.cxx",kWarning)
      << "Bad particle type for cross section shift."
      << endl;
    inlStateP = e_undefined_init_state;
    inlStateN = e_undefined_init_state;
  }
  interaction interP(e_mu,e_Fe56,e_cc,inlStateP);
  interaction interN(e_mu,e_Fe56,e_cc,inlStateP);
  Double_t defaultXSec = defaultWrapper.xsec(energy,&interP,0);
  defaultXSec += defaultWrapper.xsec(energy,&interN,0);

  neugen_cuts cuts;
  cuts.setOneProcess(e_qel);
  Double_t qelXSec = defaultWrapper.xsec(energy,&interP,&cuts);
  qelXSec += defaultWrapper.xsec(energy,&interN,&cuts);

  if (!defaultXSec){return 1.0;}
  return (defaultXSec + (this->ShiftAsValue() - 1.0)*qelXSec)/defaultXSec;
}

//____________________________________________________________________72
const Double_t NuSystematic
::ResXSecShiftScale(const Double_t energy,
                    const NuParticle::NuParticleType_t particle)
  const 
{
  neugen_config defaultConfig;
  neugen_wrapper defaultWrapper(&defaultConfig);
  init_state_t inlStateP;
  init_state_t inlStateN;
  if (NuParticle::kNuMu == particle){
    inlStateP = e_vp;
    inlStateN = e_vn;
  }
  else if (NuParticle::kNuMuBar == particle){
    inlStateP = e_vbp;
    inlStateN = e_vbn;
  }
  else{
    MSG("NuSystematic.cxx",kWarning)
      << "Bad particle type for cross section shift."
      << endl;
    inlStateP = e_undefined_init_state;
    inlStateN = e_undefined_init_state;
  }
  interaction interP(e_mu,e_Fe56,e_cc,inlStateP);
  interaction interN(e_mu,e_Fe56,e_cc,inlStateP);
  Double_t defaultXSec = defaultWrapper.xsec(energy,&interP,0);
  defaultXSec += defaultWrapper.xsec(energy,&interN,0);

  neugen_cuts cuts;
  cuts.setOneProcess(e_res);
  Double_t qelXSec = defaultWrapper.xsec(energy,&interP,&cuts);
  qelXSec += defaultWrapper.xsec(energy,&interN,&cuts);

  if (!defaultXSec){return 1.0;}
  return (defaultXSec + (this->ShiftAsValue() - 1.0)*qelXSec)/defaultXSec;
}

//____________________________________________________________________72
const Double_t NuSystematic
::NeugenXSecShiftScale(const Double_t energy,
                       const NuParticle::NuParticleType_t particle) const
{
  neugen_config shiftedConfig; //Defaults are MODBYRS4
  this->SetShiftedNeugenParameters(shiftedConfig);
  neugen_config defaultConfig;
  neugen_wrapper shiftedWrapper(&shiftedConfig);
  neugen_wrapper defaultWrapper(&defaultConfig);
  init_state_t inlStateP;
  init_state_t inlStateN;
  if (NuParticle::kNuMu == particle){
    inlStateP = e_vp;
    inlStateN = e_vn;
  }
  else if (NuParticle::kNuMuBar == particle){
    inlStateP = e_vbp;
    inlStateN = e_vbn;
  }
  else{
    MSG("NuSystematic.cxx",kWarning)
      << "Bad particle type for cross section shift."
      << endl;
    inlStateP = e_undefined_init_state;
    inlStateN = e_undefined_init_state;
  }
  interaction interP(e_mu,e_Fe56,e_cc,inlStateP);
  interaction interN(e_mu,e_Fe56,e_cc,inlStateP);
  Double_t defaultXSec = defaultWrapper.xsec(energy,&interP,0);
  defaultXSec += defaultWrapper.xsec(energy,&interN,0);
  Double_t shiftedXSec = shiftedWrapper.xsec(energy,&interP,0);
  shiftedXSec = shiftedWrapper.xsec(energy,&interN,0);
  if (defaultXSec){return shiftedXSec/defaultXSec;}
  else {return 0;}
}

---