NuSystematic.cxx

Go to the documentation of this file.

// $Id: NuSystematic.cxx,v 1.46 2010/01/30 18:59:33 bckhouse 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.46  2010/01/30 18:59:33  bckhouse
// Only apply ND cleaning systematic to NC selected events at the near detector.
//
// Change the one sigma value for CC background to 15% as used in the NC analysis.
//
// Revision 1.45  2010/01/30 00:33:36  ahimmel
//
// A significant upgrade to NuSystematic.  It can now handle multiple simultaneous shifts.  In order to handle this cleanly, significant changes were made to the internal structure of the class (e.g. how systematics and their 1 sigma values are stored).  However, it's previous functionality (apply a systematic shift based on a NuXMLConfig object) is unchanged and any previous code that uses this method should work as it did before.
//
// As with any major upgrade, some bugs might have crept in.  I've done some testing myself, but please let me know if you see some strange behavior.
//
// Here are the major new features:
//
// It is now possible to set multiple simultaneously active systematics using these methods:
//
// SetShiftsAsValues()
// SetShiftsAsSigmas()
//
// Both take an STL map with either TStrings or NuSyst::NuSystematic_t as the keys and doubles as the values.  The difference is in how the values are interpreted: SetShiftsAsValues assumes the numbers given are shifts as they might be specified in a NuXMLConfig, SetShiftsAsSigmas assumes the numbers are given in terms of sigmas.
//
// NuSystematic can also produce randomized systematic values using:
//
// Randomize()
//
// This takes any systematics that are "on" -- i.e. their current shift is not 0 and draws a random Gaussian for that systematic.
//
// There is also now a function:
//
// SetSigmas()
//
// That takes an STL map like the SetShifts functions, but now it interprets the value as the new shift that will correspond to 1 sigma.
//
// Revision 1.43  2009/09/04 15:00:06  evansj
// Fixing a bug in the resolution systematic.
//
// Revision 1.42  2009/08/28 20:48:15  gmieg
// Add '#include <cassert>'
//
// Revision 1.41  2009/08/27 14:20:40  evansj
// Added energy resolution systematics. Can shift the event energy,
// shower energy, amd track range and curvature energies independently,
// either in both detectors or just the near.
//
// Revision 1.40  2009/07/06 21:20:29  ahimmel
//
//
// Updating the DP systematic code.  FC fitting now does not consider muons part of the decay pipe.  This required adding a variable to NuFCEvent so the event library files will need to be reproduced.  Also, adjusted NuSystematic so that it used the correct codes for muons (+/-13 instead of 5/6) -- SNTPs use a different set of particle codes than the flux files.
//
// Revision 1.39  2009/07/02 14:28:33  evansj
// Don't shift muon progeny in the decay pipe systematic.
//
// Revision 1.38  2009/06/24 10:36:06  evansj
// The decay pipe systematic now only shifts events produced in the decay
// pipe (it was previously shifting those produced in the chase as well).
//
// Revision 1.37  2008/10/24 19:18:30  ahimmel
//
//
// NuFCEvent:
// Added detector to the event.  It is ok if the trees are generated with a previous version, NuFCExperimentFactory adds in the detector information when loading the events.
//
// NuFCExperimentFactory:
// Can now handle near and far detector systematics.
//
// NuMMHelperCPT:
// Updated the cross-section systematic formula.  It is not used right now, but best it be correct if we decide to later.
//
// NuMMRunFC, NuMMRunTransition, NuMMRunNuBar:
// Added a reduced constructor that does not take the far detector spectra.  This is for use in generating fake experiments with NuFCExperimentFactory.
//
// NuSystematic:
// Added a NuMuBarXSecSum that applies all appropriate systematic shifts for NuBars (both combined and ratio) using the parameterized formula that gives the sum in quadrature of all the systematic effects.
//
// Revision 1.36  2008/09/22 00:01:52  evansj
// A systematic that shifts all backgrounds in the nubar sample (CC & NC)
// by a scale factor (recommended value: 50%).
//
// Revision 1.35  2008/09/08 22:33:36  gmieg
// Change fabs(event.inu) to TMath::Abs(event.inu) to settle ambiguity.
//
// Revision 1.34  2008/08/27 12:01:42  evansj
// An uncertainty on the tau QEL & Res cross sections (simultaneaously,
// as a scale factor), the size configured from the xml files
// (recommended 15%).
//
// Revision 1.33  2008/08/19 09:44:39  evansj
// Bug fix for the target hole systematic (the hole was always being
// given a non-zero width).
//
// Revision 1.32  2008/08/18 15:11:47  evansj
// Adding a 'target hole' systematic which chops out any events in a
// segment of the target of length defined (in cm) in the xml files. The
// hole starts at z=20 cm, and is defined in terms of ppvz (which isn't
// quite correct, but my DSTs don't have tvz in them at the moment). The
// hole start of 20 cm is also a guess.
//
// Revision 1.31  2008/08/18 10:24:05  evansj
// Updating the cross section systematics so they shift all three types
// of neutrino (and antineutrino) together: not just the numus.
//
// Revision 1.30  2008/08/11 19:04:09  ahimmel
//
//
// Set up so now there is a static method in NuSystemtic called FluxSyst that returns true for flux-related systematic errors and false for detector-related ones.  This is important for the transitions analysis where the transitioned events have the detector-related systematics applied but not the flux errors.
//
// Revision 1.29  2008/06/02 11:30:38  evans
// An overall track energy scale and a near detector track energy from
// curvature scale.
//
// 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 <map>
#include <cassert>
#include <cmath>
#include "TMath.h"
#include "TString.h"
#include "TRandom3.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 "Registry/Registry.h"

#include "NtupleUtils/NuCut.h"
#include "NtupleUtils/NuEvent.h"
//#include "NtupleUtils/NuFluctuator.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuXMLConfig.h"


ClassImp(NuSystematic)

CVSID("$Id: NuSystematic.cxx,v 1.46 2010/01/30 18:59:33 bckhouse Exp $");


// Static Variables //

Bool_t NuSystematic::firstMCReweight = true;
Int_t NuSystematic::kMinusPlus = 1;
Int_t NuSystematic::kAsIs = 2;
Int_t NuSystematic::kSigma = 3;



// Constructors & Destructors //

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

  this->ConfigureDefaultSystematics();
  this->ConfigureNeugenDefaults();
  this->InitializeSystematics();
  this->EverythingOff();
  
  fCCSelector = 0;
  fNCSelector = 0;
  fNuBarSelector = 0;
  fRockSelector = 0;
  
  fRandom = new TRandom3(0);
  //fFluctuator = new NuFluctuator();
}

//____________________________________________________________________72
NuSystematic::NuSystematic(const NuXMLConfig& xmlConfig)
{
  this->ConfigureDefaultSystematics();
  this->ConfigureNeugenDefaults();
  this->InitializeSystematics();
  this->EverythingOff();
  
  fCCSelector = 0;
  fNCSelector = 0;
  fNuBarSelector = 0;
  fRockSelector = 0;
  
  fRandom = new TRandom3(0);
  //fFluctuator = new NuFluctuator();

  this->ReadXML(xmlConfig);
}

//____________________________________________________________________72
NuSystematic::~NuSystematic()
{
  //if (fFluctuator)    {delete fFluctuator;    fFluctuator = 0;}
  if (fRandom)        {delete fRandom;        fRandom = 0;}
  if (fCCSelector)    {delete fCCSelector;    fCCSelector = 0;}
  if (fNCSelector)    {delete fNCSelector;    fNCSelector = 0;}
  if (fNuBarSelector) {delete fNuBarSelector; fNuBarSelector = 0;}
  if (fRockSelector)  {delete fRockSelector;  fRockSelector = 0;}
}



// Configuration & Initialization //

//____________________________________________________________________72
void NuSystematic::ConfigureDefaultSystematics()
{
  oneSigma[NuSyst::kShowerEnergyOffset] = 100*Munits::MeV;
  oneSigma[NuSyst::kShowerEnergyScale] = 1.10;
  oneSigma[NuSyst::kShowerEnergyFunction] = 1.0;
  oneSigma[NuSyst::kShowerEnergyScaleFar] = 1.023;
  oneSigma[NuSyst::kShowerEnergyScaleNear] = 1.031;
  oneSigma[NuSyst::kShowerEnergyScaleRelative] = 1.039;
  
  oneSigma[NuSyst::kTrackEnergyCurvatureBoth] = 1.03;
  oneSigma[NuSyst::kTrackEnergyCurvatureFar] = 1.03;
  oneSigma[NuSyst::kTrackEnergyCurvatureNear] = 1.03;
  oneSigma[NuSyst::kTrackEnergyRange] = 1.02;
  oneSigma[NuSyst::kTrackEnergyScale] = 1.02;
  oneSigma[NuSyst::kTrackEnergyOverall] = -999;
  
  oneSigma[NuSyst::kBFieldBoth] = 2.0;
  oneSigma[NuSyst::kBFieldNear] = 2.0;
  oneSigma[NuSyst::kBFieldFar] = 2.0;  
  oneSigma[NuSyst::kAlignment] = -999;  
  oneSigma[NuSyst::kBeam] = 1.0;
  
  oneSigma[NuSyst::kCombinedXSecOverall] = 1.035;
  oneSigma[NuSyst::kCombinedXSecCCMA] = 1.15;
  oneSigma[NuSyst::kCombinedXSecMaRes] = 1.15;
  oneSigma[NuSyst::kCombinedXSecMaQE] = 1.15;
  oneSigma[NuSyst::kCombinedXSecDISMultip2] = 0.1;
  oneSigma[NuSyst::kCombinedXSecDISMultip3] = 0.2;
  
  oneSigma[NuSyst::kNuMuBarXSecSum] = 1.0;
  oneSigma[NuSyst::kNuMuBarXSecOverall] = 1.04;
  oneSigma[NuSyst::kNuMuBarXSecCCMA] = 1.08;
  oneSigma[NuSyst::kNuMuBarXSecDISMultip2] = 0.2;  
  oneSigma[NuSyst::kNuMuBarXSecQEL] = 1.08;
  oneSigma[NuSyst::kNuMuBarXSecRes] = 1.08;
  
  oneSigma[NuSyst::kNormalisationBoth] = 1.04;
  oneSigma[NuSyst::kNormalisationNear] = 1.04;
  oneSigma[NuSyst::kNormalisationFar] = 1.04;
  oneSigma[NuSyst::kNCBackground] = 1.50;
  oneSigma[NuSyst::kCCBackground] = 1.15;  // Value used in the NC analysis
  oneSigma[NuSyst::kNDCleaning] = 1.0;
  oneSigma[NuSyst::kAllBackgroundsScaleBoth] = 1.5;//50% scale
  oneSigma[NuSyst::kScraping] = 1.37;
  
  oneSigma[NuSyst::kJitterVDPID] = -999;
  oneSigma[NuSyst::kJitter]= 0.01;
  oneSigma[NuSyst::kDPID] = 0.03;  
  oneSigma[NuSyst::kTFProb] = 0.005;//Complete guess.
  oneSigma[NuSyst::kTargetHole] = 4.0;//cm: length of the hole
  oneSigma[NuSyst::kTauQELRes] = 1.15;//15% scale
  
  oneSigma[NuSyst::kEnergyResolutionEventBoth] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionEventNear] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionShowerBoth] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionShowerNear] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionTrackRangeBoth] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionTrackRangeNear] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionTrackCurveBoth] = 1.1;
  oneSigma[NuSyst::kEnergyResolutionTrackCurveNear] = 1.1;
}


//____________________________________________________________________72
void NuSystematic::InitializeSystematics()
{
  systNames[NuSyst::kNominal] = "Nominal";
  
  systNames[NuSyst::kShowerEnergyOffset] = "ShowerEnergyOffset";
  systNames[NuSyst::kShowerEnergyScale] = "ShowerEnergyScaleBoth"; 
  systNames[NuSyst::kShowerEnergyFunction] = "ShowerEnergyScaleFunctionBoth"; 
  systNames[NuSyst::kShowerEnergyScaleNear] = "ShowerEnergyScaleNear";
  systNames[NuSyst::kShowerEnergyScaleFar] = "ShowerEnergyScaleFar"; 
  
  systNames[NuSyst::kShowerEnergyScaleRelative] = "ShowerEnergyScaleRelative"; 
  systNames[NuSyst::kTrackEnergyCurvatureBoth] = "TrackEnergyCurvatureBoth";
  systNames[NuSyst::kTrackEnergyCurvatureFar] = "TrackEnergyCurvatureFar"; 
  systNames[NuSyst::kTrackEnergyCurvatureNear] = "TrackEnergyCurvatureNear"; 
  systNames[NuSyst::kTrackEnergyRange] = "TrackEnergyRange";
  systNames[NuSyst::kTrackEnergyOverall] = "TrackEnergyOverall"; 
  systNames[NuSyst::kTrackEnergyScale] = "TrackEnergyScale"; 
  
  systNames[NuSyst::kBFieldBoth] = "BFieldBoth"; 
  systNames[NuSyst::kBFieldNear] = "BFieldNear";
  systNames[NuSyst::kBFieldFar] = "BFieldFar";   
  systNames[NuSyst::kAlignment] = "Alignment";
  systNames[NuSyst::kBeam] = "Flux";

  systNames[NuSyst::kCombinedXSecCCMA] = "CombinedXSecCCMA"; 
  systNames[NuSyst::kCombinedXSecMaRes] = "CombinedXSecMaRes";
  systNames[NuSyst::kCombinedXSecMaQE] = "CombinedXSecMaQE"; 
  systNames[NuSyst::kCombinedXSecOverall] = "CombinedXSecOverall"; 
  systNames[NuSyst::kCombinedXSecDISMultip2] = "CombinedXSecDISMultip2";
  systNames[NuSyst::kCombinedXSecDISMultip3] = "CombinedXSecDISMultip3"; 

  systNames[NuSyst::kNuMuBarXSecCCMA] = "NuMuBarXSecCCMA"; 
  systNames[NuSyst::kNuMuBarXSecQEL] = "NuMuBarXSecQEL";
  systNames[NuSyst::kNuMuBarXSecRes] = "NuMuBarXSecRes"; 
  systNames[NuSyst::kNuMuBarXSecSum] = "NuMuBarXSecSum"; 
  systNames[NuSyst::kNuMuBarXSecOverall] = "NuMuBarXSecOverall";
  systNames[NuSyst::kNuMuBarXSecDISMultip2] = "NuMuBarXSecDISMultip2"; 
  
  systNames[NuSyst::kNormalisationBoth] = "NormalisationBoth"; 
  systNames[NuSyst::kNormalisationNear] = "NormalisationNear";
  systNames[NuSyst::kNormalisationFar] = "NormalisationFar"; 
  systNames[NuSyst::kNCBackground] = "NCBackground"; 
  systNames[NuSyst::kCCBackground] = "CCBackground";   
  systNames[NuSyst::kNDCleaning] = "NDCleaning";   
  systNames[NuSyst::kAllBackgroundsScaleBoth] = "AllBackgroundsScaleBoth"; 
  systNames[NuSyst::kScraping] = "DecayPipe"; 

  systNames[NuSyst::kJitterVDPID] = "JitterVDPID"; 
  systNames[NuSyst::kJitter] = "Jitter";
  systNames[NuSyst::kDPID] = "DPID"; 
  systNames[NuSyst::kTFProb] = "TFProb"; 
  systNames[NuSyst::kTargetHole] = "TargetHole";
  systNames[NuSyst::kTauQELRes] = "TauQELRes"; 
  
  systNames[NuSyst::kEnergyResolutionEventBoth] = "EnergyResolutionEventBoth";
  systNames[NuSyst::kEnergyResolutionEventNear] = "EnergyResolutionEventNear"; 
  systNames[NuSyst::kEnergyResolutionShowerBoth] = "EnergyResolutionShowerBoth"; 
  systNames[NuSyst::kEnergyResolutionShowerNear] = "EnergyResolutionShowerNear";
  systNames[NuSyst::kEnergyResolutionTrackRangeBoth] = "EnergyResolutionTrackRangeBoth"; 
  systNames[NuSyst::kEnergyResolutionTrackRangeNear] = "EnergyResolutionTrackRangeNear"; 
  systNames[NuSyst::kEnergyResolutionTrackCurveBoth] = "EnergyResolutionTrackCurveBoth";
  systNames[NuSyst::kEnergyResolutionTrackCurveNear] = "EnergyResolutionTrackCurveNear"; 


  systMode[NuSyst::kNominal] = kSigma;
  
  systMode[NuSyst::kShowerEnergyOffset] = kAsIs;
  systMode[NuSyst::kShowerEnergyScale] = kMinusPlus; 
  systMode[NuSyst::kShowerEnergyFunction] = kSigma; 
  systMode[NuSyst::kShowerEnergyScaleNear] = kMinusPlus;
  systMode[NuSyst::kShowerEnergyScaleFar] = kMinusPlus; 
  
  systMode[NuSyst::kShowerEnergyScaleRelative] = kMinusPlus; 
  systMode[NuSyst::kTrackEnergyCurvatureBoth] = kMinusPlus;
  systMode[NuSyst::kTrackEnergyCurvatureFar] = kMinusPlus; 
  systMode[NuSyst::kTrackEnergyCurvatureNear] = kMinusPlus; 
  systMode[NuSyst::kTrackEnergyRange] = kMinusPlus;
  systMode[NuSyst::kTrackEnergyOverall] = kSigma; 
  systMode[NuSyst::kTrackEnergyScale] = kMinusPlus; 
  
  systMode[NuSyst::kBFieldBoth] = kMinusPlus; 
  systMode[NuSyst::kBFieldNear] = kMinusPlus;
  systMode[NuSyst::kBFieldFar] = kMinusPlus;   
  systMode[NuSyst::kAlignment] = kMinusPlus;
  systMode[NuSyst::kBeam] = kSigma;
  
  systMode[NuSyst::kCombinedXSecCCMA] = kMinusPlus; 
  systMode[NuSyst::kCombinedXSecMaRes] = kMinusPlus;
  systMode[NuSyst::kCombinedXSecMaQE] = kMinusPlus; 
  systMode[NuSyst::kCombinedXSecOverall] = kMinusPlus; 
  systMode[NuSyst::kCombinedXSecDISMultip2] = kAsIs;
  systMode[NuSyst::kCombinedXSecDISMultip3] = kAsIs; 
  
  systMode[NuSyst::kNuMuBarXSecCCMA] = kMinusPlus; 
  systMode[NuSyst::kNuMuBarXSecQEL] = kMinusPlus;
  systMode[NuSyst::kNuMuBarXSecRes] = kMinusPlus; 
  systMode[NuSyst::kNuMuBarXSecSum] = kSigma; 
  systMode[NuSyst::kNuMuBarXSecOverall] = kMinusPlus;
  systMode[NuSyst::kNuMuBarXSecDISMultip2] = kAsIs; 
  
  systMode[NuSyst::kNormalisationBoth] = kMinusPlus; 
  systMode[NuSyst::kNormalisationNear] = kMinusPlus;
  systMode[NuSyst::kNormalisationFar] = kMinusPlus; 
  systMode[NuSyst::kNCBackground] = kMinusPlus; 
  systMode[NuSyst::kCCBackground] = kMinusPlus; 
  systMode[NuSyst::kNDCleaning] = kSigma; 
  systMode[NuSyst::kAllBackgroundsScaleBoth] = kMinusPlus; 
  systMode[NuSyst::kScraping] = kMinusPlus; 
  
  systMode[NuSyst::kJitterVDPID] = kSigma; 
  systMode[NuSyst::kJitter] = kAsIs;
  systMode[NuSyst::kDPID] = kAsIs; 
  systMode[NuSyst::kTFProb] = kAsIs;
  systMode[NuSyst::kTargetHole] = kAsIs;
  systMode[NuSyst::kTauQELRes] = kMinusPlus; 
  
  systMode[NuSyst::kEnergyResolutionEventBoth] = kMinusPlus;
  systMode[NuSyst::kEnergyResolutionEventNear] = kMinusPlus; 
  systMode[NuSyst::kEnergyResolutionShowerBoth] = kMinusPlus; 
  systMode[NuSyst::kEnergyResolutionShowerNear] = kMinusPlus;
  systMode[NuSyst::kEnergyResolutionTrackRangeBoth] = kMinusPlus; 
  systMode[NuSyst::kEnergyResolutionTrackRangeNear] = kMinusPlus; 
  systMode[NuSyst::kEnergyResolutionTrackCurveBoth] = kMinusPlus;
  systMode[NuSyst::kEnergyResolutionTrackCurveNear] = kMinusPlus; 
  
}


//____________________________________________________________________72
void NuSystematic::EverythingOff()
{
  currentSigma[NuSyst::kShowerEnergyOffset] = 0;
  currentSigma[NuSyst::kShowerEnergyScale] = 0;
  currentSigma[NuSyst::kShowerEnergyFunction] = 0;
  currentSigma[NuSyst::kShowerEnergyScaleFar] = 0;
  currentSigma[NuSyst::kShowerEnergyScaleNear] = 0;
  currentSigma[NuSyst::kShowerEnergyScaleRelative] = 0;
  currentSigma[NuSyst::kTrackEnergyCurvatureBoth] = 0;
  currentSigma[NuSyst::kTrackEnergyCurvatureFar] = 0;
  currentSigma[NuSyst::kTrackEnergyCurvatureNear] = 0;
  currentSigma[NuSyst::kTrackEnergyRange] = 0;
  currentSigma[NuSyst::kTrackEnergyScale] = 0;
  currentSigma[NuSyst::kTrackEnergyOverall] = 0;
  currentSigma[NuSyst::kBFieldBoth] = 0;
  currentSigma[NuSyst::kBFieldNear] = 0;
  currentSigma[NuSyst::kBFieldFar] = 0;  
  currentSigma[NuSyst::kAlignment] = 0;  
  currentSigma[NuSyst::kBeam] = 0;
  currentSigma[NuSyst::kCombinedXSecOverall] = 0;
  currentSigma[NuSyst::kCombinedXSecCCMA] = 0;
  currentSigma[NuSyst::kCombinedXSecMaRes] = 0;
  currentSigma[NuSyst::kCombinedXSecMaQE] = 0;
  currentSigma[NuSyst::kCombinedXSecDISMultip2] = 0;
  currentSigma[NuSyst::kCombinedXSecDISMultip3] = 0;
  currentSigma[NuSyst::kNuMuBarXSecSum] = 0;
  currentSigma[NuSyst::kNuMuBarXSecOverall] = 0;
  currentSigma[NuSyst::kNuMuBarXSecCCMA] = 0;
  currentSigma[NuSyst::kNuMuBarXSecDISMultip2] = 0;  
  currentSigma[NuSyst::kNuMuBarXSecQEL] = 0;
  currentSigma[NuSyst::kNuMuBarXSecRes] = 0;
  currentSigma[NuSyst::kNormalisationBoth] = 0;
  currentSigma[NuSyst::kNormalisationNear] = 0;
  currentSigma[NuSyst::kNormalisationFar] = 0;
  currentSigma[NuSyst::kNCBackground] = 0;
  currentSigma[NuSyst::kCCBackground] = 0;
  currentSigma[NuSyst::kNDCleaning] = 0;
  currentSigma[NuSyst::kScraping] = 0;
  currentSigma[NuSyst::kJitterVDPID] = 0;
  currentSigma[NuSyst::kJitter] = 0;
  currentSigma[NuSyst::kDPID] = 0;  
  currentSigma[NuSyst::kTFProb] = 0;
  currentSigma[NuSyst::kTargetHole] = 0;
  currentSigma[NuSyst::kTauQELRes] = 0;
  currentSigma[NuSyst::kAllBackgroundsScaleBoth] = 0;
  currentSigma[NuSyst::kEnergyResolutionEventBoth] = 0;
  currentSigma[NuSyst::kEnergyResolutionEventNear] = 0;
  currentSigma[NuSyst::kEnergyResolutionShowerBoth] = 0;
  currentSigma[NuSyst::kEnergyResolutionShowerNear] = 0;
  currentSigma[NuSyst::kEnergyResolutionTrackRangeBoth] = 0;
  currentSigma[NuSyst::kEnergyResolutionTrackRangeNear] = 0;
  currentSigma[NuSyst::kEnergyResolutionTrackCurveBoth] = 0;
  currentSigma[NuSyst::kEnergyResolutionTrackCurveNear] = 0;
}


//____________________________________________________________________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
void NuSystematic::ReadXML(const NuXMLConfig& xmlConfig)
{
  NuSyst::NuSystematic_t fSystematicID = SystFromName(xmlConfig.Name());
  currentSigma[fSystematicID] = 1;
  oneSigma[fSystematicID] = xmlConfig.Shift();
  
  MSG("NuSystematic.cxx",Msg::kInfo)
  << "Reading NuXMLConfig for " 
  << xmlConfig.Name() << " at "
  << currentSigma[fSystematicID] << " sigma = " << oneSigma[fSystematicID] << endl;
}

//____________________________________________________________________72
void NuSystematic::SetShiftsAsValues(map<NuSyst::NuSystematic_t, double> input)
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  MAXMSG("NuSystematic",Msg::kInfo,5) << "Setting " << input.size() 
  << " systematic shifts." << endl;
  
  for (it = input.begin(); it != input.end(); ++it) {
    currentSigma[it->first] = ConvertValueToSigma(it->second, it->first);
  }  
}


//____________________________________________________________________72
void NuSystematic::SetShiftsAsValues(map<TString, double> input)
{
  cout << "Calling SetShifts" << endl;
  map<TString, double>::const_iterator it;
  MAXMSG("NuSystematic",Msg::kInfo,5) << "Setting " << input.size() 
  << " systematic shifts." << endl;
  
  for (it = input.begin(); it != input.end(); ++it) {
    NuSyst::NuSystematic_t fSys = SystFromName(it->first);
    currentSigma[fSys] = ConvertValueToSigma(it->second, fSys);
  }  
}

//____________________________________________________________________72
void NuSystematic::SetShiftsAsSigmas(map<NuSyst::NuSystematic_t, double> input)
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  MAXMSG("NuSystematic",Msg::kInfo,5) << "Setting " << input.size() 
  << " systematic shifts." << endl;
         
  for (it = input.begin(); it != input.end(); ++it) {
    currentSigma[it->first] = it->second;
  }  
}


//____________________________________________________________________72
void NuSystematic::SetShiftsAsSigmas(map<TString, double> input)
{
  cout << "Calling SetShifts" << endl;
  map<TString, double>::const_iterator it;
  MAXMSG("NuSystematic",Msg::kInfo,5) << "Setting " << input.size() 
  << " systematic shifts." << endl;
  
  for (it = input.begin(); it != input.end(); ++it) {
    NuSyst::NuSystematic_t fSys = SystFromName(it->first);
    currentSigma[fSys] = it->second;
  }  
}


//____________________________________________________________________72
void NuSystematic::SetSigmas(map<NuSyst::NuSystematic_t, double> input)
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  
  for (it = input.begin(); it != input.end(); ++it) {
    MSG("NuSystematic",Msg::kInfo) << "Setting " << SystNames(it->first)
    << " 1 sigma to " << it->second << endl;
    oneSigma[it->first] = it->second;
  }  
  
}


//____________________________________________________________________72
void NuSystematic::SetSigmas(map<TString, double> input)
{
  map<TString, double>::const_iterator it;
  
  for (it = input.begin(); it != input.end(); ++it) {
    MSG("NuSystematic",Msg::kInfo) << "Setting " << it->first
    << " 1 sigma to " << it->second << endl;
    NuSyst::NuSystematic_t fSys = SystFromName(it->first);
    currentSigma[fSys] = it->second;
  }  
}

//____________________________________________________________________72
void NuSystematic::PrintState(bool verbose) const
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  
  
  // Find the longest name
  Int_t lsize = 12;
  for (it = currentSigma.begin(); it != currentSigma.end(); ++it) {
    if (it->second || verbose) {
      int len = SystNames(it->first).Length();
      if (len > lsize) lsize = len;
    }
  }
  
  lsize++;
  
  // The 'total width'
  // Print a summary of this events cuts
  //MSG("NuSystematic",Msg::kInfo)
  MSG("NuSystematic",Msg::kInfo) << setw(lsize+21) << left << setfill('=') << "==== NuSystematics Summary " << endl;
  
  MSG("NuSystematic",Msg::kInfo) << setw(lsize) << left << setfill(' ') << "Systematic" 
                                 << setw(10) << right << "1 sigma" << " "
                                 << setw(10) << right << "Current" << endl;
  
  for (it = currentSigma.begin(); it != currentSigma.end(); ++it) {
    if (it->second || verbose) {
      MSG("NuSystematic",Msg::kInfo) << setw(lsize) << left << setfill(' ') << SystNames(it->first) 
                                     << setw(10) << right << OneSigma(it->first) << " "
                                     << setw(10) << right << it->second << endl;
    }
  }
  MSG("NuSystematic",Msg::kInfo) << endl;
}





// Access the Systematics Maps //

//____________________________________________________________________72
NuSyst::NuSystematic_t NuSystematic::SystFromName(TString systName) const
{
  // First transform alternate names into cannonical ones
  if (systName.Contains("Scraping",TString::kIgnoreCase))
    systName = "DecayPipe";
  else if (systName.Contains("AbsoluteEnergyCalibration",TString::kIgnoreCase))
    systName = "ShowerEnergyScaleBoth";
  else if (systName.Contains("RelativeEnergyCalibrationNear",TString::kIgnoreCase))
    systName = "ShowerEnergyScaleNear";
  else if (systName.Contains("RelativeEnergyCalibrationFar",TString::kIgnoreCase))
    systName = "ShowerEnergyScaleFar";
  else if (systName.Contains("RelativeEnergyCalibration",TString::kIgnoreCase))
    systName = "ShowerEnergyScaleRelative";
  else if (systName.Contains("Beam",TString::kIgnoreCase) ||
           systName.Contains("SKZP",TString::kIgnoreCase) )
    systName = "Flux";
  else if (systName.Contains("Normalization",TString::kIgnoreCase)) {
    systName.ToLower();
    systName.ReplaceAll("lization", "lisation");
  }

  map<NuSyst::NuSystematic_t, TString>::const_iterator it;
  
  for (it = systNames.begin(); it != systNames.end(); ++it) {
    if (it->second.Contains(systName, TString::kIgnoreCase)) {
      return it->first;
    }
  }
  
  MSG("NuSystematic",Msg::kError) << "Cannot find systematic named " << systName << ", returning kUnknown." << endl;
  return NuSyst::kUnknown;
}

//____________________________________________________________________72
// The size of a 1-sigma shift
double NuSystematic::OneSigma(NuSyst::NuSystematic_t syst) const
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  it = oneSigma.find(syst);
  if (it == oneSigma.end()) {
    MSG("NuSystematic",Msg::kError) << "Systematic # " << syst << " does not exist in oneSigma." << endl;
    assert(false);
  }
  return it->second;
}

//____________________________________________________________________72
// The number of sigmas of the current shift
double NuSystematic::CurrentSigma(NuSyst::NuSystematic_t syst) const
{
  map<NuSyst::NuSystematic_t, double>::const_iterator it;
  it = currentSigma.find(syst);
  if (it == currentSigma.end()) {
    MSG("NuSystematic",Msg::kError) << "Systematic # " << syst << " does not exist in currentSigma." << endl;
    assert(false);
  }
  return it->second;  
}

//____________________________________________________________________72
// How to convert between # of sigma and a shift
const Int_t NuSystematic::SystMode(NuSyst::NuSystematic_t syst) const
{
  map<NuSyst::NuSystematic_t, Int_t>::const_iterator it;
  it = systMode.find(syst);
  if (it == systMode.end()) {
    MSG("NuSystematic",Msg::kError) << "Systematic # " << syst << " does not exist in systMode." << endl;
    assert(false);
  }
  return it->second;  
}

//____________________________________________________________________72
// The cannonical names of the systematics
const TString NuSystematic::SystNames(NuSyst::NuSystematic_t syst) const
{
  map<NuSyst::NuSystematic_t, TString>::const_iterator it;
  it = systNames.find(syst);
  if (it == systNames.end()) {
    MSG("NuSystematic",Msg::kError) << "Systematic # " << syst << " does not exist in systNames." << endl;
    assert(false);
  }
  return it->second;  
}

//____________________________________________________________________72
bool NuSystematic::FluxSyst(const NuXMLConfig& xmlConfig)
{
  TString systName = xmlConfig.Name();
  if (systName.Contains("Scraping",TString::kIgnoreCase) ||
      systName.Contains("DecayPipe",TString::kIgnoreCase)){
    return true;
  }
  if (systName.Contains("Beam",TString::kIgnoreCase) ||
      systName.Contains("SKZP",TString::kIgnoreCase) ||
      systName.Contains("Flux",TString::kIgnoreCase)){
    return true;
  }
  if (systName.Contains("TargetHole",TString::kIgnoreCase)){
    return true;
  }
  
  return false;
}



// Set the various selectors //

//____________________________________________________________________72
void NuSystematic::SetCCSelector(NuCut *input) 
{ 
  fCCSelector = input; 
  MSG("NuSystematic",Msg::kInfo) << "Setting the CC selector" << endl;
  fCCSelector->PrintSummary();
}

//____________________________________________________________________72
void NuSystematic::SetNCSelector(NuCut *input) 
{ 
  fNCSelector = input; 
  MSG("NuSystematic",Msg::kInfo) << "Setting the NC selector" << endl;
  fNCSelector->PrintSummary();
}

//____________________________________________________________________72
void NuSystematic::SetNuBarSelector(NuCut *input) 
{ 
  fNuBarSelector = input; 
  MSG("NuSystematic",Msg::kInfo) << "Setting the NuBar selector" << endl;
  fNuBarSelector->PrintSummary();
}

//____________________________________________________________________72
void NuSystematic::SetRockSelector(NuCut *input) 
{ 
  fRockSelector = input; 
  MSG("NuSystematic",Msg::kInfo) << "Setting the Rock selector" << endl;
  fRockSelector->PrintSummary();
}





// Determine Shifts and Apply Them //

//____________________________________________________________________72
void NuSystematic::Shift(NuEvent& event) const
{
  if (CurrentSigma(NuSyst::kScraping))   this->ScrapingShift(event);
  if (CurrentSigma(NuSyst::kTargetHole)) this->TargetHoleShift(event);
  if (CurrentSigma(NuSyst::kShowerEnergyOffset)) this->ShowerEnergyOffset(event);
  if (CurrentSigma(NuSyst::kShowerEnergyFunction)) this->ShowerEnergyFunction(event);
  if (CurrentSigma(NuSyst::kShowerEnergyScale) ||
      CurrentSigma(NuSyst::kShowerEnergyScaleNear) ||
      CurrentSigma(NuSyst::kShowerEnergyScaleFar) ||
      CurrentSigma(NuSyst::kShowerEnergyScaleRelative) ) this->ShowerEnergyScale(event);
  if (CurrentSigma(NuSyst::kTrackEnergyRange) ||
      CurrentSigma(NuSyst::kTrackEnergyScale) ||
      CurrentSigma(NuSyst::kTrackEnergyCurvatureBoth) ||
      CurrentSigma(NuSyst::kTrackEnergyCurvatureFar) ) this->TrackEnergyScale(event);
  if (CurrentSigma(NuSyst::kTrackEnergyOverall)) this->TrackEnergyOverall(event);
  if (CurrentSigma(NuSyst::kBeam)) this->BeamShift(event);
  if (CurrentSigma(NuSyst::kBFieldBoth) ||
      CurrentSigma(NuSyst::kBFieldNear) ||
      CurrentSigma(NuSyst::kBFieldFar) ) this->BFieldShift(event);
  if (CurrentSigma(NuSyst::kAllBackgroundsScaleBoth)) this->AllBackgroundsScaleBothShift(event);
  if (CurrentSigma(NuSyst::kNCBackground)) this->NCBackgroundShift(event);
  if (CurrentSigma(NuSyst::kCCBackground)) this->CCBackgroundShift(event);
  if (CurrentSigma(NuSyst::kNDCleaning)) this->NDCleaningShift(event);
  if (CurrentSigma(NuSyst::kNormalisationBoth) ||
      CurrentSigma(NuSyst::kNormalisationNear) ||
      CurrentSigma(NuSyst::kNormalisationFar)) this->NormalisationShift(event);
  if (CurrentSigma(NuSyst::kNuMuBarXSecSum)) this->NuMuBarSumXSecShift(event);
  if (CurrentSigma(NuSyst::kCombinedXSecOverall) ||
      CurrentSigma(NuSyst::kNuMuBarXSecOverall)) this->OverallXSecShift(event);
  if (CurrentSigma(NuSyst::kNuMuBarXSecQEL)) this->NuMuBarQELXSecShift(event);
  if (CurrentSigma(NuSyst::kNuMuBarXSecRes)) this->NuMuBarResXSecShift(event);
  if (CurrentSigma(NuSyst::kCombinedXSecCCMA) ||
      CurrentSigma(NuSyst::kCombinedXSecMaRes) ||
      CurrentSigma(NuSyst::kCombinedXSecMaQE) ||
      CurrentSigma(NuSyst::kCombinedXSecDISMultip2) ||
      CurrentSigma(NuSyst::kCombinedXSecDISMultip3) ||
      CurrentSigma(NuSyst::kNuMuBarXSecCCMA) ||
      CurrentSigma(NuSyst::kNuMuBarXSecDISMultip2)) this->NeugenXSecShift(event);
  if (CurrentSigma(NuSyst::kJitterVDPID)) this->JitterVDPIDShift(event);
  if (CurrentSigma(NuSyst::kTauQELRes)) this->TauQELResShift(event);
  if (CurrentSigma(NuSyst::kEnergyResolutionEventBoth) ||
      CurrentSigma(NuSyst::kEnergyResolutionEventNear)) this->EnergyResolutionEvent(event);
  if (CurrentSigma(NuSyst::kEnergyResolutionShowerBoth) ||
      CurrentSigma(NuSyst::kEnergyResolutionShowerNear)) this->EnergyResolutionShower(event);
  if (CurrentSigma(NuSyst::kEnergyResolutionTrackRangeBoth) ||
      CurrentSigma(NuSyst::kEnergyResolutionTrackRangeNear)) this->EnergyResolutionTrackRange(event);
  if (CurrentSigma(NuSyst::kEnergyResolutionTrackCurveBoth) ||
      CurrentSigma(NuSyst::kEnergyResolutionTrackCurveNear)) this->EnergyResolutionTrackCurve(event);
}

//____________________________________________________________________72
const Float_t NuSystematic::ConvertSigmaToValue(Float_t sigma, NuSyst::NuSystematic_t fSys) const
{  
  if (sigma == -999) sigma = CurrentSigma(fSys);
  
  if (SystMode(fSys) == kMinusPlus) {
    return sigma*(OneSigma(fSys) - 1.0) + 1.0;
  }
  else if (SystMode(fSys) == kAsIs) {
    return sigma*OneSigma(fSys);
  }
  else if (SystMode(fSys) == kSigma) {
    return sigma;
  }
  else {
    MSG("NuSystematic",Msg::kError) << "Didn't recognize syst mode " << SystMode(fSys) << " for systematic " 
    << NameFromSyst(fSys) << "(" << fSys << ")" << endl;
    return 0;
  }
}

//____________________________________________________________________72
const Float_t NuSystematic::ConvertValueToSigma(Float_t shift, NuSyst::NuSystematic_t fSys) const
{
  if (SystMode(fSys) == kMinusPlus) {
    return (shift-1.0)/(OneSigma(fSys)-1.0);
  }
  else if (SystMode(fSys) == kAsIs) {
    return shift/OneSigma(fSys);
  }
  else if (SystMode(fSys) == kSigma) {
    return shift;
  }
  else {
    MSG("NuSystematic",Msg::kError) << "Didn't recognize syst mode " << SystMode(fSys) << " for systematic " 
    << NameFromSyst(fSys) << "(" << fSys << ")" << endl;
    return 0;
  }
}

//____________________________________________________________________72
void NuSystematic::Randomize()
{
  map<NuSyst::NuSystematic_t, double>::iterator it;
  
  MAXMSG("NuSystematic",Msg::kInfo, 1) 
  << "Setting random gaussian shifts for all non-zero systematics" << endl;
  
  for (it = currentSigma.begin(); it != currentSigma.end(); ++it) {
    if (it->second) {
      it->second = fRandom->Gaus(0, 1);
      // Prevent a systematic from being turned off 
      // persistently by a random draw
      if (!it->second) it->second = 1e-8; 
    }
  }
}



// Apply systematic shifts //

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

//____________________________________________________________________72
void NuSystematic::JitterVDPIDShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track jitter & DPID shift = " << ShiftAsSigma(NuSyst::kJitterVDPID) << "s" << endl;
  event.jitter += ConvertSigmaToValue(ShiftAsSigma(NuSyst::kJitterVDPID), NuSyst::kJitter);
  event.dpID += ConvertSigmaToValue(ShiftAsSigma(NuSyst::kJitterVDPID), NuSyst::kDPID);
  return;
}

//____________________________________________________________________72
void NuSystematic::TFProbShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track fit probability shift = " << ShiftAsSigma(NuSyst::kTFProb) << "s" << endl;
  event.prob += ShiftAsValue(NuSyst::kTFProb);
  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;
  
  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, event);

  
  //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::NuMuBarSumXSecShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing numubar summed cross section shift = "
    << ShiftAsSigma(NuSyst::kNuMuBarXSecSum) << "s" << endl;
    if (1 != event.iaction) return;
    if (event.inu != -14) return;
    
    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 = event.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 + ShiftAsValue(NuSyst::kNuMuBarXSecSum) * shift;
    
    event.rw *= shift;
    return;
}

//____________________________________________________________________72
void NuSystematic::OverallXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing overall cross section shift combined = " << ShiftAsSigma(NuSyst::kCombinedXSecOverall)
  << "s, numubar = " << ShiftAsSigma(NuSyst::kNuMuBarXSecOverall) << "s" << endl;
  if (1 != event.iaction) return;
  
  event.rw *= ShiftAsValue(NuSyst::kCombinedXSecOverall);
  if (event.inu < 0) 
    event.rw *= ShiftAsValue(NuSyst::kNuMuBarXSecOverall);
}

//____________________________________________________________________72
void NuSystematic::TauQELResShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing tau QEL+Res shift: "
    << (this->ShiftAsValue(NuSyst::kTauQELRes)-1.)*100. << "%." << endl;

  if (16 == TMath::Abs(event.inu) && 1 == event.iaction
      && (1001 == event.iresonance || 1002 == event.iresonance)){
    event.rw *= ShiftAsValue(NuSyst::kTauQELRes);
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::NuMuBarQELXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NuMuBar QEL cross section shift = "
    << ShiftAsSigma(NuSyst::kNuMuBarXSecQEL) << "s" << endl;

  if (1 != event.iaction) return;
  //  if (-14 != event.inu) return;
  if (event.inu > 0) return;
  if (1001 != event.iresonance) return;
  event.rw *= ShiftAsValue(NuSyst::kNuMuBarXSecQEL);
  return;
}

//____________________________________________________________________72
void NuSystematic::NuMuBarResXSecShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NuMuBar Res cross section shift"
    << ShiftAsSigma(NuSyst::kNuMuBarXSecRes) << "s" << endl;

  if (1 != event.iaction) return;
  //  if (-14 != event.inu) return;
  if (event.inu > 0) return;
  if (1002 != event.iresonance) return;
  event.rw *= ShiftAsValue(NuSyst::kNuMuBarXSecRes);
  return;
}

//____________________________________________________________________72
void NuSystematic::AllBackgroundsScaleBothShift(NuEvent& event) const
{
  if (!fNuBarSelector) {
    MSG("NuSystematic",Msg::kError) << "Cannot apply AllBackgroundsScaleBoth without doing SetNuBarSelector()" << endl;
    assert(false);
  }
  
  fNuBarSelector->MakeCuts(event);
  if (1==event.charge && fNuBarSelector->Passed()
      && (14 == event.inu || 0 == event.iaction)) {

    MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing AllBackgroundsScaleBoth shift ="
    << ShiftAsSigma(NuSyst::kAllBackgroundsScaleBoth) << "s" << endl;
    
    event.rw *= ShiftAsValue(NuSyst::kAllBackgroundsScaleBoth);
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::BFieldShift(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing BField shift both = " << ShiftAsSigma(NuSyst::kBFieldBoth) 
    << "s, near = " << ShiftAsSigma(NuSyst::kBFieldNear) 
    << "s, far = " << ShiftAsSigma(NuSyst::kBFieldFar) << endl;
  if (1==event.charge && 14 == event.inu && 1 == event.iaction){
    event.rw *= ShiftAsValue(NuSyst::kBFieldBoth);
    if (Detector::kNear==event.detector)
      event.rw *= ShiftAsValue(NuSyst::kBFieldNear);
    if (Detector::kFar==event.detector)
      event.rw *= ShiftAsValue(NuSyst::kBFieldFar);
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::NCBackgroundShift(NuEvent& event) const
{
  if (!fCCSelector) {
    MSG("NuSystematic",Msg::kError) << "Cannot apply NCBackground without doing SetCCSelector()" << endl;
    assert(false);
  }
  
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing NC background shift " << ShiftAsSigma(NuSyst::kNCBackground) << "s" << endl;
  
  fCCSelector->MakeCuts(event);
  if (0 == event.iaction && fCCSelector->Passed()) 
    event.rw *= ShiftAsValue(NuSyst::kNCBackground);
}

//____________________________________________________________________72
void NuSystematic::CCBackgroundShift(NuEvent& event) const
{
  if (!fNCSelector) {
    MSG("NuSystematic",Msg::kError) << "Cannot apply CCBackground without doing SetNCSelector()" << endl;
    assert(false);
  }
  
  MAXMSG("NuSystematic",Msg::kInfo,1)
  << "Performing CC background shift " << ShiftAsSigma(NuSyst::kCCBackground) << "s" << endl;
  
  fNCSelector->MakeCuts(event);
  if (1 == event.iaction && fNCSelector->Passed()) 
    event.rw *= ShiftAsValue(NuSyst::kCCBackground);
}

//____________________________________________________________________72
void NuSystematic::NDCleaningShift(NuEvent& event) const
{
  if (!fNCSelector) {
    MSG("NuSystematic",Msg::kError) << "Cannot apply NDCleaning without doing SetNCSelector()" << endl;
    assert(false);
  }
  
  MAXMSG("NuSystematic",Msg::kInfo,1)
  << "Performing ND cleaning shift " << ShiftAsSigma(NuSyst::kNDCleaning) << "s" << endl;
  
  double percent = 0;
  if(event.energyNC < 0.5) percent = 15.2;
  if(event.energyNC >= 0.5 && event.energyNC < 1.0) percent = 2.9;
  if(event.energyNC >= 1.0 && event.energyNC < 1.5) percent = 0.4;
  
  fNCSelector->MakeCuts(event);
  if(fNCSelector->Passed() && event.detector == Detector::kNear)
    event.rw *= 1.+(percent/100.)*ShiftAsValue(NuSyst::kNDCleaning);
}

//____________________________________________________________________72
void NuSystematic::NormalisationShift(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing Normalisation shift both = "  << ShiftAsSigma(NuSyst::kNormalisationBoth)
      << "s, near = " << ShiftAsSigma(NuSyst::kNormalisationNear)
      << "s, far = " << ShiftAsSigma(NuSyst::kNormalisationFar) << endl;
  
  event.rw *= ShiftAsValue(NuSyst::kNormalisationBoth);
  if (Detector::kNear==event.detector)
    event.rw *= ShiftAsValue(NuSyst::kNormalisationNear);
  if (Detector::kFar==event.detector)
    event.rw *= ShiftAsValue(NuSyst::kNormalisationFar);
  return;
}

//____________________________________________________________________72
void NuSystematic::TargetHoleShift(NuEvent& event) const
{
  Double_t holeStart = 20.0;
  Double_t holeLength = ShiftAsValue(NuSyst::kTargetHole);

  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing target hole shift. Hole length: "
    << holeLength << " cm" << endl;

  if ((event.ppvz > holeStart) && (event.ppvz < holeStart+holeLength)){
    event.rw = 0;
  }
  return;
}

//____________________________________________________________________72
void NuSystematic::ScrapingShift(NuEvent& event) const 
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing scraping shift = " << ShiftAsSigma(NuSyst::kScraping) << "s" << 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;
//    }

    //Don't shift muons
    if (13 == event.ptype) return;
    if (-13 == event.ptype) return;
    
//    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

    if (event.ppvz < 4500)  return;
    //Not decay pipe
    
    event.rw *= ShiftAsValue(NuSyst::kScraping);
    
    return;
}

//____________________________________________________________________72
void NuSystematic::ShowerEnergyOffset(NuEvent& event) const
{
    MAXMSG("NuSystematic",Msg::kInfo,1)
      << "Performing shower energy offset shift = " << ShiftAsValue(NuSyst::kShowerEnergyOffset)/Munits::MeV << " MeV" << endl;
  event.shwEn += ShiftAsValue(NuSyst::kShowerEnergyOffset)/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 " 
    << ShiftAsSigma(NuSyst::kShowerEnergyFunction)
    << "s" << endl;

  Double_t offset = 7.0;
  Double_t scale = 4.0;
  Double_t eLife = 1.5; //GeV
  Double_t enShift = 1.0 + this->ShiftAsSigma(NuSyst::kShowerEnergyFunction)*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 with near = " << ShiftAsSigma(NuSyst::kShowerEnergyScaleNear)
      << "s, far = " << ShiftAsSigma(NuSyst::kShowerEnergyScaleFar)
      << "s, relative = " << ShiftAsSigma(NuSyst::kShowerEnergyScaleRelative)
      << "s, overall = " << ShiftAsSigma(NuSyst::kShowerEnergyScale) << "s" << endl;
  
  if (Detector::kNear==event.detector) {
    event.shwEn *= ShiftAsValue(NuSyst::kShowerEnergyScaleNear);
    event.shwEn *= ShiftAsValue(NuSyst::kShowerEnergyScaleRelative);
  }
  if (Detector::kFar==event.detector) {
    event.shwEn *= ShiftAsValue(NuSyst::kShowerEnergyScaleFar);
  }
  event.shwEn *= ShiftAsValue(NuSyst::kShowerEnergyScale);
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::TrackEnergyScale(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
  << "Performing track energy scale with range = " << ShiftAsSigma(NuSyst::kTrackEnergyRange)
  << "s, curve near = " << ShiftAsSigma(NuSyst::kTrackEnergyCurvatureNear)
  << "s, curve far = " << ShiftAsSigma(NuSyst::kTrackEnergyCurvatureFar)
  << "s, curve both = " << ShiftAsSigma(NuSyst::kTrackEnergyCurvatureBoth)
  << "s, scale = " << ShiftAsSigma(NuSyst::kShowerEnergyScale) << "s" << endl;
  
  if (event.usedRange){
    event.trkEn *= ShiftAsValue(NuSyst::kTrackEnergyRange);
  }
  else {
    event.trkEn *= ShiftAsValue(NuSyst::kTrackEnergyCurvatureBoth);
    if (Detector::kNear==event.detector)
      event.trkEn *= ShiftAsValue(NuSyst::kTrackEnergyCurvatureNear);
    if (Detector::kFar==event.detector)
      event.trkEn *= ShiftAsValue(NuSyst::kTrackEnergyCurvatureFar);
  }
  event.trkEn *= ShiftAsValue(NuSyst::kTrackEnergyScale);
  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(NuSyst::kTrackEnergyOverall) << " 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(NuSyst::kTrackEnergyOverall)*0.01*eScale;
  event.energy = event.shwEn+event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::EnergyResolutionEvent(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing event energy resolution shift both = " << ShiftAsSigma(NuSyst::kEnergyResolutionEventBoth)
    << "s, near = " << ShiftAsSigma(NuSyst::kEnergyResolutionEventNear) << "s" << endl;
  
  if (fabs(event.neuEnMC) < 1e-9) return;
  Double_t deltaE = (event.energy - event.neuEnMC) / event.neuEnMC;
  deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionEventBoth);
  if (Detector::kNear == event.detector)
    deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionEventNear);  
  event.energy = deltaE*event.neuEnMC + event.neuEnMC;
  return;
}

//____________________________________________________________________72
void NuSystematic::EnergyResolutionShower(NuEvent& event) const
{
  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing shower energy resolution shift both = " << ShiftAsSigma(NuSyst::kEnergyResolutionShowerBoth)
    << "s, near = " << ShiftAsSigma(NuSyst::kEnergyResolutionShowerNear) << "s" << endl;
  
  if (fabs(event.shwEnMC) <= 1e-9) return;
  if (event.nshw < 1) return;
  if (fabs(event.shwEn) < 1e-9) return;
  Double_t deltaE = (event.shwEn - event.shwEnMC) / event.shwEnMC;
  deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionShowerBoth);
  if (Detector::kNear == event.detector)
    deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionShowerNear);  
  event.shwEn = deltaE*event.shwEnMC + event.shwEnMC;
  event.energy = event.shwEn + event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::EnergyResolutionTrackRange(NuEvent& event) const
{
  if (!event.usedRange) return;

  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track energy range resolution shift both = " << ShiftAsSigma(NuSyst::kEnergyResolutionTrackRangeBoth)
    << "s, near = " << ShiftAsSigma(NuSyst::kEnergyResolutionTrackRangeNear) << "s" << endl;

  if (event.trkEnRange != event.trkEn){
    cout << "Stopping event, but event.trkEnRange != event.trkEn"
    << endl;
    assert(false);
  }

  if (fabs(event.trkEnMC) <= 1e-9) return;
  if (event.ntrk < 1) return;
  if (fabs(event.trkEn) < 1e-9) return;
  Double_t deltaE = (event.trkEnRange - event.trkEnMC) / event.trkEnMC;
  deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionTrackRangeBoth);
  if (Detector::kNear == event.detector)
    deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionTrackRangeNear);  
  
  event.trkEnRange = deltaE*event.trkEnMC + event.trkEnMC;

  event.trkEn = event.trkEnRange;
  event.energy = event.shwEn + event.trkEn;
  return;
}

//____________________________________________________________________72
void NuSystematic::EnergyResolutionTrackCurve(NuEvent& event) const
{
  if (!event.usedCurv) return;

  MAXMSG("NuSystematic",Msg::kInfo,1)
    << "Performing track energy curve resolution shift both = " << ShiftAsSigma(NuSyst::kEnergyResolutionTrackCurveBoth)
    << "s, near = " << ShiftAsSigma(NuSyst::kEnergyResolutionTrackCurveNear) << "s" << endl;
  
  
  if (event.trkEnCurv != event.trkEn){
    cout << "Exiting event, but event.trkEnCurv != event.trkEn"
    << endl;
    assert(false);
  }


  if (fabs(event.trkEnMC) <= 1e-9) return;
  if (event.ntrk < 1) return;
  if (fabs(event.trkEn) < 1e-9) return;
  Double_t deltaE = (event.trkEnCurv - event.trkEnMC) / event.trkEnMC;
  deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionTrackCurveBoth);
  if (Detector::kNear == event.detector)
    deltaE *= ShiftAsValue(NuSyst::kEnergyResolutionTrackCurveNear);  

  event.trkEnCurv = deltaE*event.trkEnMC + event.trkEnMC;

  event.trkEn = event.trkEnCurv;
  event.energy = event.shwEn + event.trkEn;
  return;
}





// Neugen Functions //

//____________________________________________________________________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::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 NuEvent event) const
{
  Float_t ma_qe = this->MA_QEDefault();
  Float_t ma_res = this->MA_ResDefault();
  
  ma_qe *= ShiftAsValue(NuSyst::kCombinedXSecCCMA);
  ma_res *= ShiftAsValue(NuSyst::kCombinedXSecCCMA);
  ma_qe *= ShiftAsValue(NuSyst::kCombinedXSecMaQE);
  ma_res *= ShiftAsValue(NuSyst::kCombinedXSecMaRes);
  
  if (event.inu<0) {
    ma_qe *= ShiftAsValue(NuSyst::kNuMuBarXSecCCMA);
    ma_res *= ShiftAsValue(NuSyst::kNuMuBarXSecCCMA);
  }

  registry.Set("neugen:ma_qe",ma_qe);
  registry.Set("neugen:ma_res",ma_res);

  Float_t kno112 = this->KNO112Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno112<0.0){kno112=0.0;}
  if (kno112>1.0){kno112=1.0;}
  Float_t kno122 = this->KNO122Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno122<0.0){kno122=0.0;}
  if (kno122>1.0){kno122=1.0;}
  Float_t kno132 = this->KNO132Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (event.inu<0) kno132 += this->ShiftAsValue(NuSyst::kNuMuBarXSecDISMultip2);
  if (kno132<0.0){kno132=0.0;}
  if (kno132>1.0){kno132=1.0;}
  Float_t kno142 = this->KNO142Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (event.inu<0) kno142 += this->ShiftAsValue(NuSyst::kNuMuBarXSecDISMultip2);
  if (kno142<0.0){kno142=0.0;}
  if (kno142>1.0){kno142=1.0;}
  Float_t kno212 = this->KNO212Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno212<0.0){kno212=0.0;}
  if (kno212>1.0){kno212=1.0;}
  Float_t kno222 = this->KNO222Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno222<0.0){kno222=0.0;}
  if (kno222>1.0){kno222=1.0;}
  Float_t kno232 = this->KNO232Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno232<0.0){kno232=0.0;}
  if (kno232>1.0){kno232=1.0;}
  Float_t kno242 = this->KNO242Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  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);

  Float_t kno113 = this->KNO113Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno113<0.0){kno113=0.0;}
  if (kno113>1.0){kno113=1.0;}
  Float_t kno123 = this->KNO123Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno123<0.0){kno123=0.0;}
  if (kno123>1.0){kno123=1.0;}
  Float_t kno133 = this->KNO133Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno133<0.0){kno133=0.0;}
  if (kno133>1.0){kno133=1.0;}
  Float_t kno143 = this->KNO143Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno143<0.0){kno143=0.0;}
  if (kno143>1.0){kno143=1.0;}
  Float_t kno213 = this->KNO213Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno213<0.0){kno213=0.0;}
  if (kno213>1.0){kno213=1.0;}
  Float_t kno223 = this->KNO223Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno223<0.0){kno223=0.0;}
  if (kno223>1.0){kno223=1.0;}
  Float_t kno233 = this->KNO233Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno233<0.0){kno233=0.0;}
  if (kno233>1.0){kno233=1.0;}
  Float_t kno243 = this->KNO243Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  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);
}

//____________________________________________________________________72
void NuSystematic::SetShiftedNeugenParameters(neugen_config& config, const NuEvent event) const
{
  Float_t ma_qe = this->MA_QEDefault();
  Float_t ma_res = this->MA_ResDefault();
  
  ma_qe *= ShiftAsValue(NuSyst::kCombinedXSecCCMA);
  ma_res *= ShiftAsValue(NuSyst::kCombinedXSecCCMA);
  ma_qe *= ShiftAsValue(NuSyst::kCombinedXSecMaQE);
  ma_res *= ShiftAsValue(NuSyst::kCombinedXSecMaRes);
  
  if (event.inu<0) {
    ma_qe *= ShiftAsValue(NuSyst::kNuMuBarXSecCCMA);
    ma_res *= ShiftAsValue(NuSyst::kNuMuBarXSecCCMA);
  }
  
  config.set_ma_qe(ma_qe);
  config.set_ma_res(ma_res);
  
  Float_t kno112 = this->KNO112Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno112<0.0){kno112=0.0;}
  if (kno112>1.0){kno112=1.0;}
  Float_t kno122 = this->KNO122Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno122<0.0){kno122=0.0;}
  if (kno122>1.0){kno122=1.0;}
  Float_t kno132 = this->KNO132Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (event.inu<0) kno132 += this->ShiftAsValue(NuSyst::kNuMuBarXSecDISMultip2);
  if (kno132<0.0){kno132=0.0;}
  if (kno132>1.0){kno132=1.0;}
  Float_t kno142 = this->KNO142Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (event.inu<0) kno142 += this->ShiftAsValue(NuSyst::kNuMuBarXSecDISMultip2);
  if (kno142<0.0){kno142=0.0;}
  if (kno142>1.0){kno142=1.0;}
  Float_t kno212 = this->KNO212Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno212<0.0){kno212=0.0;}
  if (kno212>1.0){kno212=1.0;}
  Float_t kno222 = this->KNO222Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno222<0.0){kno222=0.0;}
  if (kno222>1.0){kno222=1.0;}
  Float_t kno232 = this->KNO232Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  if (kno232<0.0){kno232=0.0;}
  if (kno232>1.0){kno232=1.0;}
  Float_t kno242 = this->KNO242Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip2);
  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);
  
  Float_t kno113 = this->KNO113Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno113<0.0){kno113=0.0;}
  if (kno113>1.0){kno113=1.0;}
  Float_t kno123 = this->KNO123Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno123<0.0){kno123=0.0;}
  if (kno123>1.0){kno123=1.0;}
  Float_t kno133 = this->KNO133Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno133<0.0){kno133=0.0;}
  if (kno133>1.0){kno133=1.0;}
  Float_t kno143 = this->KNO143Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno143<0.0){kno143=0.0;}
  if (kno143>1.0){kno143=1.0;}
  Float_t kno213 = this->KNO213Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno213<0.0){kno213=0.0;}
  if (kno213>1.0){kno213=1.0;}
  Float_t kno223 = this->KNO223Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno223<0.0){kno223=0.0;}
  if (kno223>1.0){kno223=1.0;}
  Float_t kno233 = this->KNO233Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  if (kno233<0.0){kno233=0.0;}
  if (kno233>1.0){kno233=1.0;}
  Float_t kno243 = this->KNO243Default() + this->ShiftAsValue(NuSyst::kCombinedXSecDISMultip3);
  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);
}




// Xsec Shift Functions //

//____________________________________________________________________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 (CurrentSigma(NuSyst::kCombinedXSecOverall)){
    return this->ShiftAsValue(NuSyst::kCombinedXSecOverall);
  }
  else if (CurrentSigma(NuSyst::kNuMuBarXSecOverall)){
    if (NuParticle::kNuMuBar == particle){
      return this->ShiftAsValue(NuSyst::kNuMuBarXSecOverall);
    }
    else {
      return 1.0;
    }
  }
  else if (CurrentSigma(NuSyst::kNuMuBarXSecCCMA) ||
           CurrentSigma(NuSyst::kNuMuBarXSecDISMultip2)){
    if (NuParticle::kNuMuBar != particle){
      return 1.0;
    }
    else {return this->NeugenXSecShiftScale(energy,particle);}
  }
  else if(CurrentSigma(NuSyst::kCombinedXSecCCMA) ||
          CurrentSigma(NuSyst::kCombinedXSecDISMultip2) ||
          CurrentSigma(NuSyst::kCombinedXSecDISMultip3)){
    if (!(NuParticle::kNuMu == particle || NuParticle::kNuMuBar == particle)){
      return 1.0;
    }
    else if(CurrentSigma(NuSyst::kNuMuBarXSecQEL)){
      return this->QELXSecShiftScale(energy,particle);
    }
    else if (CurrentSigma(NuSyst::kNuMuBarXSecRes)){
      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(NuSyst::kNuMuBarXSecQEL) - 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(NuSyst::kNuMuBarXSecRes) - 1.0)*qelXSec)/defaultXSec;
}

//____________________________________________________________________72
const Double_t NuSystematic
::NeugenXSecShiftScale(const Double_t energy,
                       const NuParticle::NuParticleType_t particle) const
{
  NuEvent event;
  if (NuParticle::kNuMuBar == particle) event.inu = -14;
  else                                  event.inu = 14;

  neugen_config shiftedConfig; //Defaults are MODBYRS4
  this->SetShiftedNeugenParameters(shiftedConfig, event);
  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;}
}

---