NuFluxHelper.cxx

Go to the documentation of this file.

// $Id: NuFluxHelper.cxx,v 1.39 2010/02/08 22:52:39 ahimmel Exp $
//
// Makes fluka-related helper files for Matrix Method extrapolation.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
// Adapted from code by Chris Smith
//
//
// $Log: NuFluxHelper.cxx,v $
// Revision 1.39  2010/02/08 22:52:39  ahimmel
//
// Use the new kDogwood1_Daikon07_v2 weights.
//
// Revision 1.38  2010/01/15 00:11:35  ahimmel
//
//
// Set reproducible random seeds based on run number of the first file.  If that number cannot be determined, default to using the root-based unique seed.
//
// Revision 1.37  2010/01/13 21:27:10  ahimmel
//
// Now set up for parallelized production of flux files.
//
// Revision 1.36  2009/12/19 19:14:57  ahimmel
//
// Fixed how filenames are determined for choosing whether you have a Flugg or Gnumi tree. Now globs rather than using the wildcarded path.
//
// Revision 1.35  2009/11/18 00:43:58  evansj
// Creating a beam matrix which converts from an ND numu+numubar flux X
// cross section to a FD nutau+nutaubar flux X cross section.
//
// Revision 1.34  2009/10/26 22:50:02  ahimmel
//
// Upgrades for using new Dogwood1_Daikon07 SKZP weights added.  Now automatically uses kDetXs for Daikon04 and earlier and kDogwood1_Daikon07 for Daikon07+
//
// Added a new unified way to print out the release type that prints: Name (hex #, decimal #).
//
// Revision 1.33  2009/10/02 22:04:33  ahimmel
//
//
// Added a way to define the particles to extrapolate that doesn't use std::vector so it can play nice with CINT.
//
// Revision 1.32  2009/09/29 19:40:27  ahimmel
//
//
// Was picking up the beam weights but not using them.
//
// Revision 1.31  2009/09/13 16:37:14  nickd
// Added missing <cassert> include to fix compile
//
// Alex' last check in does not compile without it on Mac platform
//
// Revision 1.30  2009/09/12 22:58:37  ahimmel
//
//
// The beam matrix making code is now up-to-date to handle the new Flugg fluxes.  It does intelligent things when run periods/beam types are specified or not, double checking where possible and complaining about errors.  It should make sure that only the right files are used to make the desired beam matrix.
//
// Revision 1.29  2009/08/25 15:10:05  evansj
// Updating the beam matrix implementation of NUWEIGHT to recognise the
// omega and k0short.
//
// Revision 1.28  2009/08/25 14:28:20  evansj
// The new Flugg files store data members as Double_t, whereas the old
// Gnumi files stored them as Float_t. This means all the
// SetBranchAddress commends in NuFluxChain no longer work. So
// NuFluxChain is now a virtual base class. NuFluggChain and NuGnumiChain
// inherit from it (using doubles and floats respectively). The beam
// matrix code looks at the flux file name to see which of these classes
// it should instantiate.
//
// Revision 1.27  2008/03/21 14:56:59  evans
// Adding some printout to reassure the user the reweighting is being
// done for the correct run and beam type.
//
// Revision 1.26  2008/03/21 12:59:13  evans
// NuFluxHelper is now configured from the xml file to perform SKZP and
// scraping systematics shifts.
//
// Revision 1.25  2008/03/03 17:51:41  evans
// Commenting out the line of code that makes the non MM-reweighted beam
// matrices, as these are never needed. This speeds the code up.
//
// Revision 1.24  2008/02/22 11:12:01  evans
// Updating the positions of the near and far detectors in the beam
// matrix code to agree with Robert Hatcher's Daikon_04 numbers.
//
// Revision 1.23  2008/02/20 19:15:31  evans
// Changing the ND position for the pHE beam. Adding some useful kInfo
// printout.
//
// Revision 1.22  2008/02/19 16:10:21  evans
// Changing the name of the new function in NuZBeamReweight from
// GetZBeamReweight to GetBeamWeightsOnly.
//
// Revision 1.21  2008/02/19 15:21:40  evans
// A new function in NuZBeamReweight to get beam-only SKZP weights for
// the beam matrix.
//
// Revision 1.20  2008/02/19 12:15:50  evans
// The beam matrix binning is now configured from the .xml file. Default
// is 0.5 GeV bins up to 200 GeV.
//
// Revision 1.19  2008/02/19 11:50:55  evans
// Updating the SKZP reweighting code so it now uses the NuZBeamReweight
// interface.
//
// Revision 1.18  2008/02/08 21:49:21  evans
// Updating the position of the near detector.
//
// Revision 1.17  2008/01/23 22:40:56  evans
// A new constructor taking an xml filename (as a TString) that
// configures whether or not to SKZP, and which run period to use.
//
// Revision 1.16  2008/01/22 21:04:25  hartnell
//
// Used to use:
//
// nuCuts.IsInFidVolNDCC0325Std(...)
//
// but now use:
//
// nuCuts.IsInFidVol(...)
//
// Revision 1.15  2008/01/21 01:16:13  evans
// Changing to the new CC0325 ND fiducial volume.
//
// Revision 1.14  2008/01/15 01:47:10  evans
// Bug fix: was subtracting 1 from SKZP errors.
//
// Revision 1.13  2008/01/14 07:00:37  evans
// Beam matrices can now be made with a decay-pipe scraping systematic
// shift. Also allowing the SKZP systematic to be used event if the beam
// matrix isn't bein SKZP-reweighted (not entirely sure if that's what we
// want or not).
//
// Revision 1.12  2008/01/14 06:39:22  evans
// The beam matrix can now be made systematically shifted by the SKZP
// flux errors.
//
// Revision 1.11  2008/01/08 20:33:13  evans
// Setting the pointer 'fzarko' to 0 at the start of the constructor, so
// that when the function DoSKZP(Bool_t) tries to construct it the test
// 'if (!fzarko)' confirms the pointer is not yet pointing to
// anything. This prevents a seg fault on a compiler in which pointers
// aren't, by default, set to 0.
//
// Revision 1.10  2007/12/14 23:51:27  ahimmel
//
//
// Made some modifications so that when run w/o skzp reweighting, the skzp object is not constructed (fixed some crashes when running with questionable db access that caused producing the skzp object to crash).
//
// Revision 1.9  2007/11/11 06:42:06  rhatcher
// Part of the ongoing purge of "DetectorType" in favor of "Detector"
// (the former is a synonym of the latter ...)
//
// Revision 1.8  2007/09/22 17:09:33  evans
// *** empty log message ***
//
// Revision 1.7  2007/09/22 16:55:45  evans
// Adding the ability to do SKZP reweighting for run I & run II combined.
//
// Revision 1.6  2007/08/07 16:42:19  evans
// Beam matrices can now be made for any ND fiducial volume (needn't be
// cylindrical).
//
// This involved altering the ND co-ordinate to be the real ND co-ordinate
// system (no longer centered on the beam spot), and making the ND fiducial
// volume and position of the beam spot centre more easily configurable. At
// present the old beam spot center is still hard coded in, but can easily be
// changed.
//
// Revision 1.5  2007/08/06 16:15:13  evans
// Allowing more than one type of neutrino to be extrapolated at once. Only
// CC numus and CC numubars are supported at the moment.
//
// The macro MakeMMFluxHelpers.C must now be compiled as the use of the
// NuParticle enum doesn't work in interpreted mode.
//
// Revision 1.4  2007/08/03 10:16:11  evans
// Fixing a bug in the setting of branch addresses in NuFluxChain.
//
// Fixing a bug in how one of the matrices was normalised in NuFluxHelper.
//
// Fixing the units where NuFluxChain (which uses Munits) interacts with the
// matrix method code in NuFluxHelper (which doesn't use Munits).
//
// Now making only one TRandom3 object as a member variable of NuFluxHelper so
// that the random neutrino interaction point in the near detector really is
// random.
//
// Revision 1.3  2007/08/02 15:19:49  evans
// In NuHistos: filling the data histograms (RecoEnergy_ND & _FD) for both
// data and MC, so that mock data can be used.
//
// Improving the memory management in NuFluxHelper: primarily removing the
// memory leak whereby ten new TRandom3 objects were created for each fluka
// entry and never deleted.
//
// Revision 1.2  2007/08/01 04:43:23  rhatcher
// need #include <cmath> for fabs() on gcc 4.0.1.
//
// Revision 1.1  2007/07/26 10:32:14  evans
// Adding an implementation of code to make a beam matrix for extrapolation.
//
// NuFluxHelper makes the beam matrix. It's configurable to make the matrix
// for either numus or numubars, and even to make a non-square beam matrix if
// the mood takes you.
//
// NuFluxChain is a wrapper around a fluka file to give the interface required
// by NuFluxHelper.
//
// macros/MakeMMFluxHelpers.C is an example macro of how to make the beam
// matrix.
//
//
//

#include <cmath>
#include <vector>
#include <cassert>

#include "TFile.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH1F.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TString.h"

#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuCuts.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuFluggChain.h"
#include "NtupleUtils/NuFluxChain.h"
#include "NtupleUtils/NuFluxHelper.h"
#include "NtupleUtils/NuGnumiChain.h"
#include "NtupleUtils/NuXMLConfig.h"
#include "NtupleUtils/NuZBeamReweight.h"
#include "NtupleUtils/NuUtilities.h"
#include "Conventions/SimFlag.h"

using std::string;
using std::vector;

using NuParticle::NuParticleType_t;

ClassImp(NuFluxHelper)

CVSID("$Id: NuFluxHelper.cxx,v 1.39 2010/02/08 22:52:39 ahimmel Exp $");


//____________________________________________________________________72
NuFluxHelper::NuFluxHelper(SKZPWeightCalculator::ERunPeriod rp, Bool_t _doRW)
{
  batchRunning = false;
  fzarko = 0;
  DoSKZP(_doRW);
  frunPeriod = rp;
  //fbeamType = BeamType::kL010z185i;
  fbeamType = BeamType::kUnknown; // Now pulled from file by default
  freweightVersion = SKZPWeightCalculator::kUnknown; // Set when determining file version

  fdoBeamShift = false;
  fdoScrapingShift = false;
  fbeamShiftAsSigma = 0.0;
  fscrapingScaleFactor = 0.0;

  fbinningScheme = NuBinningScheme::kUnknown;
  const NuUtilities utils;  
  const vector<Double_t> trueBins = utils.TrueBins(fbinningScheme);
  fNNDTrueEBins = trueBins.size()-1;
  fNDTrueEBinEdges=new Double_t[fNNDTrueEBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = trueBins.begin();
         itBin != trueBins.end();
         ++itBin, ++i){
      fNDTrueEBinEdges[i] = *itBin;
    }
  }
  fNFDTrueEBins = fNNDTrueEBins;
  fFDTrueEBinEdges = fNDTrueEBinEdges;

  foutFile = "";
  fxSecFile = "";

  fxSecGraphNuMuCC = 0;
  fxSecGraphNuMuBarCC = 0;
  fxSecGraphNuTauCC = 0;
  fxSecGraphNuTauBarCC = 0;

  fFDVsNDMatrix = 0;        
  fFDVsNDMatrixRW = 0;                              
  fFDVsNDMatrixXSec = 0;    
  fFDVsNDMatrixXSecRW = 0;  
  fFDVsNDMatrixTauXSecRW = 0;
                        
  fTrueEnergyNuFlux_ND = 0;
  fTrueEnergyNuFluxRW_ND = 0;
  fTrueEnergyNuFlux_FD = 0;
  fTrueEnergyNuFluxRW_FD = 0;
  fTrueEnergyCCFlux_ND = 0;
  fTrueEnergyCCFluxRW_ND = 0;
  fTrueEnergyCCFlux_FD = 0;
  fTrueEnergyCCFluxRW_FD = 0;

  frandom3 = new TRandom3();
}

//____________________________________________________________________72
NuFluxHelper::NuFluxHelper(const TString& xmlFile)
{
  cout << "Beginning to construct NuFluxHelper(" << xmlFile << ")" << endl;
  fzarko = 0;
  batchRunning = false;

  NuXMLConfig xmlConfig(xmlFile);
  if (-1==xmlConfig.RunPeriod()){
    frunPeriod = SKZPWeightCalculator::kNone;
    this->DoSKZP(false);
    cout << "NOT doing SKZP." << endl;
  }
  else if (0==xmlConfig.RunPeriod()){
    frunPeriod = SKZPWeightCalculator::kNone;
    this->DoSKZP(true);
  }
  else if (1==xmlConfig.RunPeriod()){
    frunPeriod = SKZPWeightCalculator::kRunI;
    this->DoSKZP(true);
  }
  else if (2==xmlConfig.RunPeriod()){
    frunPeriod = SKZPWeightCalculator::kRunII;
    this->DoSKZP(true);
  }
  else {
    frunPeriod = SKZPWeightCalculator::kNone;
    this->DoSKZP(false);
  }
  cout << "Picked a run period (" << xmlConfig.RunPeriod() << ")" << endl;
  
  //fbeamType = BeamType::kL010z185i;
  fbeamType = BeamType::kUnknown; // Now pulled from file by default
  freweightVersion = SKZPWeightCalculator::kUnknown; // Set when determining file version
  
  if (xmlConfig.Name().Contains("Beam",TString::kIgnoreCase) ||
      xmlConfig.Name().Contains("SKZP",TString::kIgnoreCase) ||
      xmlConfig.Name().Contains("Flux",TString::kIgnoreCase)){
    this->SystematicBeamShift(true,xmlConfig.Shift());
  }
  else{
    this->SystematicBeamShift(false, 0.0);
  }
  if (xmlConfig.Name().Contains("Scraping",TString::kIgnoreCase) ||
      xmlConfig.Name().Contains("DecayPipe",TString::kIgnoreCase)){
    this->SystematicScrapingShift(true,xmlConfig.Shift());
  }
  else{
    this->SystematicScrapingShift(false, 0.0);
  }
  cout << "Read other things from the xml file" << endl;
  
  fbinningScheme = static_cast<NuBinningScheme::NuBinningScheme_t>(xmlConfig.BinningScheme());
  const NuUtilities utils;
  
  const vector<Double_t> trueBins = utils.TrueBins(fbinningScheme);
  fNNDTrueEBins = trueBins.size()-1;
  fNDTrueEBinEdges=new Double_t[fNNDTrueEBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = trueBins.begin();
         itBin != trueBins.end();
         ++itBin, ++i){
      fNDTrueEBinEdges[i] = *itBin;
    }
  }
  fNFDTrueEBins = fNNDTrueEBins;
  fFDTrueEBinEdges = fNDTrueEBinEdges;
  
  cout << "Constructed bins" << endl;
  
  foutFile = "";
  fxSecFile = "";

  fxSecGraphNuMuCC = 0;
  fxSecGraphNuMuBarCC = 0;
  fxSecGraphNuTauCC = 0;
  fxSecGraphNuTauBarCC = 0;

  fFDVsNDMatrix = 0;        
  fFDVsNDMatrixRW = 0;                              
  fFDVsNDMatrixXSec = 0;    
  fFDVsNDMatrixXSecRW = 0;  
  fFDVsNDMatrixTauXSecRW = 0;
                        
  fTrueEnergyNuFlux_ND = 0;
  fTrueEnergyNuFluxRW_ND = 0;
  fTrueEnergyNuFlux_FD = 0;
  fTrueEnergyNuFluxRW_FD = 0;
  fTrueEnergyCCFlux_ND = 0;
  fTrueEnergyCCFluxRW_ND = 0;
  fTrueEnergyCCFlux_FD = 0;
  fTrueEnergyCCFluxRW_FD = 0;

  frandom3 = new TRandom3();
  
  cout << "Initialized objects" << endl;
}

//____________________________________________________________________72
void NuFluxHelper::DoSKZP(const Bool_t doSKZP)
{
  fdoSKZP=doSKZP;
  if (fdoSKZP && !fzarko) {
    fzarko = new NuZBeamReweight();
  }
}

//____________________________________________________________________72
NuFluxHelper::~NuFluxHelper()
{
  if (fFDVsNDMatrix)
    {delete fFDVsNDMatrix; fFDVsNDMatrix = 0;}
  if (fFDVsNDMatrixRW)
    {delete fFDVsNDMatrixRW; fFDVsNDMatrixRW = 0;}
  if (fFDVsNDMatrixXSec)
    {delete fFDVsNDMatrixXSec; fFDVsNDMatrixXSec = 0;}
  if (fFDVsNDMatrixXSecRW)
    {delete fFDVsNDMatrixXSecRW; fFDVsNDMatrixXSecRW = 0;}
  if (fFDVsNDMatrixTauXSecRW)
    {delete fFDVsNDMatrixTauXSecRW; fFDVsNDMatrixTauXSecRW = 0;}
  
  if (fTrueEnergyNuFlux_ND)
    {delete fTrueEnergyNuFlux_ND; fTrueEnergyNuFlux_ND = 0;}
  if (fTrueEnergyNuFluxRW_ND)
    {delete fTrueEnergyNuFluxRW_ND; fTrueEnergyNuFluxRW_ND = 0;}
  if (fTrueEnergyNuFlux_FD)
    {delete fTrueEnergyNuFlux_FD; fTrueEnergyNuFlux_FD = 0;}
  if (fTrueEnergyNuFluxRW_FD)
    {delete fTrueEnergyNuFluxRW_FD; fTrueEnergyNuFluxRW_FD = 0;}
  if (fTrueEnergyCCFlux_ND)
    {delete fTrueEnergyCCFlux_ND; fTrueEnergyCCFlux_ND = 0;}
  if (fTrueEnergyCCFluxRW_ND)
    {delete fTrueEnergyCCFluxRW_ND; fTrueEnergyCCFluxRW_ND = 0;}
  if (fTrueEnergyCCFlux_FD)
    {delete fTrueEnergyCCFlux_FD; fTrueEnergyCCFlux_FD = 0;}
  if (fTrueEnergyCCFluxRW_FD)
    {delete fTrueEnergyCCFluxRW_FD; fTrueEnergyCCFluxRW_FD = 0;}

  if (fFDTrueEBinEdges) {delete[] fFDTrueEBinEdges; fFDTrueEBinEdges=0;}
  if (fNDTrueEBinEdges) {delete[] fNDTrueEBinEdges; fNDTrueEBinEdges=0;}

  if (fzarko) {delete fzarko; fzarko = 0;}

  if (fxSecGraphNuMuCC) {delete fxSecGraphNuMuCC; fxSecGraphNuMuCC = 0;}
  if (fxSecGraphNuMuBarCC) {delete fxSecGraphNuMuBarCC; fxSecGraphNuMuBarCC = 0;}
  if (fxSecGraphNuTauCC) {delete fxSecGraphNuTauCC; fxSecGraphNuTauCC = 0;}
  if (fxSecGraphNuTauBarCC) {delete fxSecGraphNuTauBarCC; fxSecGraphNuTauBarCC = 0;}

  if (frandom3) {delete frandom3; frandom3 = 0;}
}

//____________________________________________________________________72
void NuFluxHelper::FDTrueEBins(const vector<Double_t> fdtruebins)
{
  Int_t size = fdtruebins.size();
  fNFDTrueEBins = size - 1;
  if (fFDTrueEBinEdges) {delete[] fFDTrueEBinEdges; fFDTrueEBinEdges=0;}
  fFDTrueEBinEdges = new Double_t[size];
  Int_t i = 0;
  for (vector<Double_t>::const_iterator it = fdtruebins.begin();
       it != fdtruebins.end();
       ++it, ++i){
    fFDTrueEBinEdges[i] = *it;
  }
}

//____________________________________________________________________72
void NuFluxHelper::NDTrueEBins(const vector<Double_t> ndtruebins)
{
  Int_t size = ndtruebins.size();
  fNNDTrueEBins = size - 1;
  if (fNDTrueEBinEdges) {delete[] fNDTrueEBinEdges; fNDTrueEBinEdges=0;}
  fNDTrueEBinEdges = new Double_t[size];
  Int_t i = 0;
  for (vector<Double_t>::const_iterator it = ndtruebins.begin();
       it != ndtruebins.end();
       ++it, ++i){
    fNDTrueEBinEdges[i] = *it;
  }
}

//____________________________________________________________________72
const Bool_t NuFluxHelper::CreateHistos()
{
  if (fFDVsNDMatrix)
    {delete fFDVsNDMatrix; fFDVsNDMatrix = 0;}
  if (fFDVsNDMatrixRW)
    {delete fFDVsNDMatrixRW; fFDVsNDMatrixRW = 0;}
  if (fFDVsNDMatrixXSec)
    {delete fFDVsNDMatrixXSec; fFDVsNDMatrixXSec = 0;}
  if (fFDVsNDMatrixXSecRW)
    {delete fFDVsNDMatrixXSecRW; fFDVsNDMatrixXSecRW = 0;}
  if (fFDVsNDMatrixTauXSecRW)
    {delete fFDVsNDMatrixTauXSecRW; fFDVsNDMatrixTauXSecRW = 0;}
  
  if (fTrueEnergyNuFlux_ND)
    {delete fTrueEnergyNuFlux_ND; fTrueEnergyNuFlux_ND = 0;}
  if (fTrueEnergyNuFluxRW_ND)
    {delete fTrueEnergyNuFluxRW_ND; fTrueEnergyNuFluxRW_ND = 0;}
  if (fTrueEnergyNuFlux_FD)
    {delete fTrueEnergyNuFlux_FD; fTrueEnergyNuFlux_FD = 0;}
  if (fTrueEnergyNuFluxRW_FD)
    {delete fTrueEnergyNuFluxRW_FD; fTrueEnergyNuFluxRW_FD = 0;}
  if (fTrueEnergyCCFlux_ND)
    {delete fTrueEnergyCCFlux_ND; fTrueEnergyCCFlux_ND = 0;}
  if (fTrueEnergyCCFluxRW_ND)
    {delete fTrueEnergyCCFluxRW_ND; fTrueEnergyCCFluxRW_ND = 0;}
  if (fTrueEnergyCCFlux_FD)
    {delete fTrueEnergyCCFlux_FD; fTrueEnergyCCFlux_FD = 0;}
  if (fTrueEnergyCCFluxRW_FD)
    {delete fTrueEnergyCCFluxRW_FD; fTrueEnergyCCFluxRW_FD = 0;}
  
  if (!fNNDTrueEBins){
    MSG("NuFluxHelper",Msg::kError)
      <<"Number of ND true energy bins not set " <<endl;
    return false;
  }
  if (!fNFDTrueEBins){
    MSG("NuFluxHelper",Msg::kError)
      <<"Number of FD true energy bins not set " <<endl;
    return false;
  }
  if (!fNDTrueEBinEdges){
    MSG("NuFluxHelper",Msg::kError)
      <<"ND true energy bin edges not set " <<endl;
    return false;
  }
  if (!fFDTrueEBinEdges){
    MSG("NuFluxHelper",Msg::kError)
      <<"FD true energy bin edges not set " <<endl;
    return false;
  }
  
  fFDVsNDMatrix = new 
    TH2D("FDVsNDMatrix",
         "Number of FD Vs ND Events with True Energy",
         fNNDTrueEBins,fNDTrueEBinEdges,
         fNFDTrueEBins,fFDTrueEBinEdges);
  fFDVsNDMatrix->Sumw2();
  fFDVsNDMatrixRW = new 
    TH2D("FDVsNDMatrixRW",
         "Number of FD Vs ND Events with True Energy (with Near Reweight)",
         fNNDTrueEBins,fNDTrueEBinEdges,
         fNFDTrueEBins,fFDTrueEBinEdges);
  fFDVsNDMatrixRW->Sumw2();
  fFDVsNDMatrixXSec = new 
    TH2D("FDVsNDMatrixXSec",
         "Number of FD Vs ND Events with True Energy (with XSec)",
         fNNDTrueEBins,fNDTrueEBinEdges,
         fNFDTrueEBins,fFDTrueEBinEdges);
  fFDVsNDMatrixXSec->Sumw2();
  fFDVsNDMatrixXSecRW = new 
    TH2D("FDVsNDMatrixXSecRW",
         "Number of FD Vs ND Events with True Energy (with XSec + Near Reweight)",
         fNNDTrueEBins,fNDTrueEBinEdges,
         fNFDTrueEBins,fFDTrueEBinEdges);
  fFDVsNDMatrixXSecRW->Sumw2();
  fFDVsNDMatrixTauXSecRW = new
    TH2D("FDVsNDMatrixTauXSecRW",
         "Number of FD Vs ND Events with True Energy (with FD Tau XSec + Near Reweight)",
         fNNDTrueEBins,fNDTrueEBinEdges,
         fNFDTrueEBins,fFDTrueEBinEdges);
  
  fTrueEnergyNuFlux_ND = new 
    TH1D("TrueEnergyNuFlux_ND",
         "Neutrino Flux with True Energy (NearDet)",
         fNNDTrueEBins,fNDTrueEBinEdges);
  fTrueEnergyNuFlux_ND->Sumw2();
  fTrueEnergyNuFluxRW_ND = new 
    TH1D("TrueEnergyNuFluxRW_ND",
         "Neutrino Flux with True Energy (NearDet with Reweighting)",
         fNNDTrueEBins,fNDTrueEBinEdges);
  fTrueEnergyNuFlux_FD = new 
    TH1D("TrueEnergyNuFlux_FD",
         "Neutrino Flux with True Energy (FarDet)",
         fNFDTrueEBins,fFDTrueEBinEdges);
  fTrueEnergyNuFlux_FD->Sumw2();
  fTrueEnergyNuFluxRW_FD = new 
    TH1D("TrueEnergyNuFluxRW_FD",
         "Neutrino Flux with True Energy (FarDet with Reweighting)",
         fNFDTrueEBins,fFDTrueEBinEdges);
  fTrueEnergyNuFluxRW_FD->Sumw2();
  fTrueEnergyCCFlux_ND = new 
    TH1D("TrueEnergyCCFlux_ND",
         "NuMu CC Flux with True Energy (NearDet)",
         fNNDTrueEBins,fNDTrueEBinEdges);
  fTrueEnergyCCFlux_ND->Sumw2();
  fTrueEnergyCCFluxRW_ND = new
    TH1D("TrueEnergyCCFluxRW_ND",
         "NuMu CC Flux with True Energy with (NearDet with Reweighting)",
         fNNDTrueEBins,fNDTrueEBinEdges);
  fTrueEnergyCCFluxRW_ND->Sumw2();  
  fTrueEnergyCCFlux_FD = new
    TH1D("TrueEnergyCCFlux_FD",
         "NuMu CC Flux with True Energy (FarDet)",
         fNFDTrueEBins,fFDTrueEBinEdges);
  fTrueEnergyCCFlux_FD->Sumw2();
  fTrueEnergyCCFluxRW_FD = new
    TH1D("TrueEnergyCCFluxRW_FD",
         "NuMu CC Flux with True Energy (FarDet with Reweighing)",
         fNFDTrueEBins,fFDTrueEBinEdges);
  fTrueEnergyCCFluxRW_FD->Sumw2();
  
  return true;
}

//____________________________________________________________________72
void NuFluxHelper::CrossSectionFile(const char * xSecFile)
{
  fxSecFile = xSecFile;

  if (fxSecGraphNuMuCC){
    delete fxSecGraphNuMuCC;
    fxSecGraphNuMuCC = 0;
  }
  if (fxSecGraphNuMuBarCC){
    delete fxSecGraphNuMuBarCC;
    fxSecGraphNuMuBarCC = 0;
  }
  if (fxSecGraphNuTauCC){
    delete fxSecGraphNuTauCC;
    fxSecGraphNuTauCC = 0;
  }
  if (fxSecGraphNuTauBarCC){
    delete fxSecGraphNuTauBarCC;
    fxSecGraphNuTauBarCC = 0;
  }
  
  //Get the cross-section information from the histogram in the file.
  //Put it into an array.
  TFile *mikeFile=new TFile(fxSecFile.c_str(),"READ");

  //Put the numu CC cross sections into an array
  TH1F* XSec_NuMuCC = (TH1F*) mikeFile->Get("h_numu_cc_tot");
  XSec_NuMuCC->SetDirectory(0);
  Float_t* x = new Float_t[XSec_NuMuCC->GetNbinsX()];
  Float_t* y = new Float_t[XSec_NuMuCC->GetNbinsX()];
  for(int i=0;i<XSec_NuMuCC->GetNbinsX();i++) {
    x[i] = XSec_NuMuCC->GetBinCenter(i+1);
    y[i] = XSec_NuMuCC->GetBinContent(i+1);
  }
  //Turn the array into a graph.
  fxSecGraphNuMuCC = new TGraph(XSec_NuMuCC->GetNbinsX(),x,y);

  //Clean up ready for the next set of cross sections
  if (x) {delete x; x=0;}
  if (y) {delete y; y=0;}

  //Put the numubar CC cross sections into an array
  TH1F* XSec_NuMuBarCC = (TH1F*) mikeFile->Get("h_numubar_cc_tot");
  XSec_NuMuBarCC->SetDirectory(0);
  x = new Float_t[XSec_NuMuBarCC->GetNbinsX()];
  y = new Float_t[XSec_NuMuBarCC->GetNbinsX()];
  for(int i=0;i<XSec_NuMuBarCC->GetNbinsX();i++) {
    x[i] = XSec_NuMuBarCC->GetBinCenter(i+1);
    y[i] = XSec_NuMuBarCC->GetBinContent(i+1);
  }
  //Turn the array into a graph.
  fxSecGraphNuMuBarCC = new TGraph(XSec_NuMuBarCC->GetNbinsX(),x,y);

  //Clean up ready for the next set of cross sections:
  if (x) {delete x; x=0;}
  if (y) {delete y; y=0;}

  //Put the nutau CC cross sections into an array
  TH1F* XSec_NuTauCC = (TH1F*) mikeFile->Get("h_nutau_cc_tot");
  XSec_NuTauCC->SetDirectory(0);
  x = new Float_t[XSec_NuTauCC->GetNbinsX()];
  y = new Float_t[XSec_NuTauCC->GetNbinsX()];
  for(int i=0;i<XSec_NuTauCC->GetNbinsX();i++) {
    x[i] = XSec_NuTauCC->GetBinCenter(i+1);
    y[i] = XSec_NuTauCC->GetBinContent(i+1);
  }
  //Turn the array into a graph.
  fxSecGraphNuTauCC = new TGraph(XSec_NuTauCC->GetNbinsX(),x,y);

  //Clean up ready for the next set of cross sections
  if (x) {delete x; x=0;}
  if (y) {delete y; y=0;}

  //Put the nutaubar CC cross sections into an array
  TH1F* XSec_NuTauBarCC = (TH1F*) mikeFile->Get("h_nutaubar_cc_tot");
  XSec_NuTauBarCC->SetDirectory(0);
  x = new Float_t[XSec_NuTauBarCC->GetNbinsX()];
  y = new Float_t[XSec_NuTauBarCC->GetNbinsX()];
  for(int i=0;i<XSec_NuTauBarCC->GetNbinsX();i++) {
    x[i] = XSec_NuTauBarCC->GetBinCenter(i+1);
    y[i] = XSec_NuTauBarCC->GetBinContent(i+1);
  }
  //Turn the array into a graph.
  fxSecGraphNuTauBarCC = new TGraph(XSec_NuTauBarCC->GetNbinsX(),x,y);

  //Clean up memory:
  if (x) {delete x; x=0;}
  if (y) {delete y; y=0;}

  mikeFile->Close();
  if (mikeFile) {delete mikeFile; mikeFile = 0;}
}

//____________________________________________________________________72
void NuFluxHelper::
ParticlesToExtrapolate(int flavour)
{
  //Configure the particles to extrapolate
  vector<NuParticle::NuParticleType_t> particles;
  if (0==flavour){
    particles.push_back(NuParticle::kNuMu);
    particles.push_back(NuParticle::kNuMuBar);
  }
  if (1==flavour){
    particles.push_back(NuParticle::kNuMu);
  }
  if (2==flavour){
    particles.push_back(NuParticle::kNuMuBar);
  }

  this->ParticlesToExtrapolate(particles);
}


//____________________________________________________________________72
void NuFluxHelper::
ParticleToExtrapolate(const NuParticleType_t particle)
{
  fparticlesToExtrapolate.clear();
  fparticlesToExtrapolate.push_back(particle);
}

//____________________________________________________________________72
void NuFluxHelper::
ParticlesToExtrapolate(const vector<NuParticleType_t> particleList)
{
  fparticlesToExtrapolate.clear();
  fparticlesToExtrapolate = particleList;
}

//____________________________________________________________________72
const Bool_t NuFluxHelper::
IsParticleToExtrapolate(const NuParticleType_t particle) const
{
  for (vector<NuParticleType_t>::const_iterator it
         = fparticlesToExtrapolate.begin();
       it != fparticlesToExtrapolate.end();
       ++it){
    if (particle == *it){return true;}
  }
  return false;
}

//____________________________________________________________________72
void NuFluxHelper::MakeHelperHistos(const char * fileList)
{
  const NuFluxChain* fluxChain = 0;
  vector<TString> vFiles = NuUtilities::GetListOfFilesInDir(fileList);
  TString tmp = vFiles[0];
  if (tmp.Contains("flugg")){
    MSG("NuFluxHelper",Msg::kInfo)  << "USING FLUGG CHAIN" << endl;
    fluxChain = new NuFluggChain(fileList);
    // Reset the reweight version only for Flugg
    freweightVersion = SKZPWeightCalculator::kDogwood1_Daikon07_v2;
  }
  else{
    MSG("NuFluxHelper",Msg::kInfo) << "USING GNUMI CHAIN" << endl;
    fluxChain = new NuGnumiChain(fileList);
    freweightVersion = SKZPWeightCalculator::kDetXs;
  }

  TString num = tmp(tmp.Last('_')+1, 3);
  if (!num.IsDigit()) {
    MSG("NuFluxHelper",Msg::kInfo) << "Cannot determine run number for file: " << tmp
    << " (got " << num << ") so random seed will not be reproducible." << endl;
    frandom3->SetSeed(0);
  }
  else {
    MSG("NuFluxHelper",Msg::kInfo) << "Setting random seed to " << num << "." << endl;
    frandom3->SetSeed(num.Atoi());
  }
    
  
  
  // Double check the beam type
  // If our beam type is not yet known, fill it with type from fluxChain
  // If fluxChain type is not known, trust our beam type
  // If we have both types, check that they match
  if (fluxChain->GetBeamType() != BeamType::kUnknown) {
    if (BeamType() == BeamType::kUnknown) BeamType(fluxChain->GetBeamType());
    if (BeamType() != fluxChain->GetBeamType()) {
      MSG("NuFluxHelper",Msg::kError) << "Beam types do not match.  "
      << "Type from fluxFiles: " << BeamType::AsString(fluxChain->GetBeamType())
      << ", Type from FluxHelper: " << BeamType::AsString(BeamType()) << endl;
      return;
    }
  }
  
  

  // Double check the Run Period
  // If our run period is not yet known, fill it with type from fluxChain
  // If fluxChain run period is not known, trust our run period
  // If we have both types, check that they match
  if (fluxChain->GetRunPeriod() != SKZPWeightCalculator::kNone) {
    if (RunPeriod() == SKZPWeightCalculator::kNone) RunPeriod(fluxChain->GetRunPeriod());
    if (RunPeriod() != fluxChain->GetRunPeriod()) {
      MSG("NuFluxHelper",Msg::kError) << "Run Periods do not match.  "
      << "Period from fluxFiles: " << static_cast<int>(fluxChain->GetRunPeriod())
      << ", Period from FluxHelper: " << static_cast<int>(RunPeriod()) << endl;
      return;
    }
  }
  
  MSG("NuFluxHelper",Msg::kInfo) << "Make flux helpers with " << BeamType::AsString(BeamType())
                                 << " and run period " << static_cast<int>(RunPeriod()) << endl;
  
  //Make the histograms if possible:
  Bool_t histosReady = this->CreateHistos();
  if (!histosReady){return;}

  if (!fxSecGraphNuMuCC){
    MSG("NuFluxHelper",Msg::kError)
      <<"No numu CC cross-section graph supplied." <<endl;
    return;
  }
  if (!fxSecGraphNuMuBarCC){
    MSG("NuFluxHelper",Msg::kError)
      <<"No numubar CC cross-section graph supplied." <<endl;
    return;
  }
  if (!fxSecGraphNuTauCC){
    MSG("NuFluxHelper",Msg::kError)
      <<"No nutau CC cross-section graph supplied." <<endl;
    return;
  }
  if (!fxSecGraphNuTauBarCC){
    MSG("NuFluxHelper",Msg::kError)
      <<"No nutaubar CC cross-section graph supplied." <<endl;
    return;
  }

  //Loop over the entries in the fluka files.
  Int_t numFlukaEntries = fluxChain->GetEntries();
  for(int flukaEntry=0;flukaEntry<numFlukaEntries;++flukaEntry){
    if(flukaEntry%10000==0){
      MSG("NuFluxHelper",Msg::kInfo)
        <<"Processing fluka entry " << flukaEntry
        << " of " << numFlukaEntries <<endl;
    }

    //Is this the correct particle?
    if (!this->IsParticleToExtrapolate(fluxChain->NuType(flukaEntry))){
      continue;
    }
    /*
    this->FillNonMMHelpers(fluxChain,
                           flukaEntry);
    */
    this->FillMMHelpers(fluxChain,
                        flukaEntry);

  } //End loop over flukaChain

  if (!batchRunning) {
    this->NormaliseHistos();
  }
  this->WriteHistos();

  if (fluxChain){delete fluxChain; fluxChain=0;}

  return;
}


//____________________________________________________________________72
void NuFluxHelper::ConcatenateHelpers(const char * fileList)
{
  vector<TString> vFiles = NuUtilities::GetListOfFilesInDir(fileList);
  if (vFiles.size() <= 0) {
    MSG("NuFluxHelper",Msg::kInfo) << "Cannot find any files in " << fileList << endl;
  }
  
  MSG("NuFluxHelper",Msg::kInfo) << "Concatenating " << vFiles.size()
  << " flux helper files." << endl;
  
  TFile *f;
  
  MSG("NuFluxHelper",Msg::kInfo) << "Reading from " << vFiles[0] << endl;
  f = new TFile(vFiles[0]);
  if (!f->IsOpen()) {
    MSG("NuFluxHelper",Msg::kError) << "File not opened.  Bailing" << endl;
    assert(false);
  }
  f->ls();
  
  fFDVsNDMatrix = (TH2D*)f->Get("FDVsNDMatrix"); 
  fFDVsNDMatrixRW = (TH2D*)f->Get("FDVsNDMatrixRW"); 
  fFDVsNDMatrixXSec = (TH2D*)f->Get("FDVsNDMatrixXSec"); 
  fFDVsNDMatrixXSecRW = (TH2D*)f->Get("FDVsNDMatrixXSecRW"); 
  fFDVsNDMatrixTauXSecRW = (TH2D*)f->Get("FDVsNDMatrixTauXSecRW"); 
  
  fTrueEnergyNuFlux_ND = (TH1D*)f->Get("TrueEnergyNuFlux_ND"); 
  fTrueEnergyNuFluxRW_ND = (TH1D*)f->Get("TrueEnergyNuFluxRW_ND"); 
  fTrueEnergyNuFlux_FD = (TH1D*)f->Get("TrueEnergyNuFlux_FD"); 
  fTrueEnergyNuFluxRW_FD = (TH1D*)f->Get("TrueEnergyNuFluxRW_FD"); 
  fTrueEnergyCCFlux_ND = (TH1D*)f->Get("TrueEnergyCCFlux_ND"); 
  fTrueEnergyCCFluxRW_ND = (TH1D*)f->Get("TrueEnergyCCFluxRW_ND"); 
  fTrueEnergyCCFlux_FD = (TH1D*)f->Get("TrueEnergyCCFlux_FD"); 
  fTrueEnergyCCFluxRW_FD = (TH1D*)f->Get("TrueEnergyCCFluxRW_FD"); 

  fFDVsNDMatrix->SetDirectory(0);
  fFDVsNDMatrixRW->SetDirectory(0);
  fFDVsNDMatrixXSec->SetDirectory(0);
  fFDVsNDMatrixXSecRW->SetDirectory(0);
  fFDVsNDMatrixTauXSecRW->SetDirectory(0);
  
  fTrueEnergyNuFlux_ND->SetDirectory(0);
  fTrueEnergyNuFluxRW_ND->SetDirectory(0);
  fTrueEnergyNuFlux_FD->SetDirectory(0);
  fTrueEnergyNuFluxRW_FD->SetDirectory(0);
  fTrueEnergyCCFlux_ND->SetDirectory(0);
  fTrueEnergyCCFluxRW_ND->SetDirectory(0);
  fTrueEnergyCCFlux_FD->SetDirectory(0);
  fTrueEnergyCCFluxRW_FD->SetDirectory(0);
  f->Close();
  
  
  for (unsigned int i = 1; i < vFiles.size(); i++) {
    MSG("NuFluxHelper",Msg::kInfo) << "Reading from " << vFiles[i] << endl;
    f = new TFile(vFiles[i]);
    if (!f->IsOpen()) {
      MSG("NuFluxHelper",Msg::kError) << "File not opened.  Bailing" << endl;
      assert(false);
    }
    
    TH2D* tmp_fFDVsNDMatrix = (TH2D*)f->Get("FDVsNDMatrix");
    TH2D* tmp_fFDVsNDMatrixRW = (TH2D*)f->Get("FDVsNDMatrixRW");
    TH2D* tmp_fFDVsNDMatrixXSec = (TH2D*)f->Get("FDVsNDMatrixXSec");
    TH2D* tmp_fFDVsNDMatrixXSecRW = (TH2D*)f->Get("FDVsNDMatrixXSecRW");
    TH2D* tmp_fFDVsNDMatrixTauXSecRW = (TH2D*)f->Get("FDVsNDMatrixTauXSecRW");
    
    TH1D* tmp_fTrueEnergyNuFlux_ND = (TH1D*)f->Get("TrueEnergyNuFlux_ND");
    TH1D* tmp_fTrueEnergyNuFluxRW_ND = (TH1D*)f->Get("TrueEnergyNuFluxRW_ND");
    TH1D* tmp_fTrueEnergyNuFlux_FD = (TH1D*)f->Get("TrueEnergyNuFlux_FD");
    TH1D* tmp_fTrueEnergyNuFluxRW_FD = (TH1D*)f->Get("TrueEnergyNuFluxRW_FD");
    TH1D* tmp_fTrueEnergyCCFlux_ND = (TH1D*)f->Get("TrueEnergyCCFlux_ND");
    TH1D* tmp_fTrueEnergyCCFluxRW_ND = (TH1D*)f->Get("TrueEnergyCCFluxRW_ND");
    TH1D* tmp_fTrueEnergyCCFlux_FD = (TH1D*)f->Get("TrueEnergyCCFlux_FD");
    TH1D* tmp_fTrueEnergyCCFluxRW_FD = (TH1D*)f->Get("TrueEnergyCCFluxRW_FD");

    fFDVsNDMatrix->Add(tmp_fFDVsNDMatrix);
    fFDVsNDMatrixRW->Add(tmp_fFDVsNDMatrixRW);
    fFDVsNDMatrixXSec->Add(tmp_fFDVsNDMatrixXSec);
    fFDVsNDMatrixXSecRW->Add(tmp_fFDVsNDMatrixXSecRW);
    fFDVsNDMatrixTauXSecRW->Add(tmp_fFDVsNDMatrixTauXSecRW);
    
    fTrueEnergyNuFlux_ND->Add(tmp_fTrueEnergyNuFlux_ND);
    fTrueEnergyNuFluxRW_ND->Add(tmp_fTrueEnergyNuFluxRW_ND);
    fTrueEnergyNuFlux_FD->Add(tmp_fTrueEnergyNuFlux_FD);
    fTrueEnergyNuFluxRW_FD->Add(tmp_fTrueEnergyNuFluxRW_FD);
    fTrueEnergyCCFlux_ND->Add(tmp_fTrueEnergyCCFlux_ND);
    fTrueEnergyCCFluxRW_ND->Add(tmp_fTrueEnergyCCFluxRW_ND);
    fTrueEnergyCCFlux_FD->Add(tmp_fTrueEnergyCCFlux_FD);
    fTrueEnergyCCFluxRW_FD->Add(tmp_fTrueEnergyCCFluxRW_FD);

    f->Close();
  }
  
  this->NormaliseHistos();
  this->WriteHistos();
  
  return;
}

//____________________________________________________________________72
void NuFluxHelper::WriteHistos() const
{
  MSG("NuFluxHelper",Msg::kInfo)
    <<"Writing helper histograms to " << foutFile << endl;
  
  TFile *savefile = new TFile(foutFile.c_str(),"RECREATE");
  savefile->cd();
  fFDVsNDMatrix->Write();
  fFDVsNDMatrixRW->Write();
  fFDVsNDMatrixXSec->Write();
  fFDVsNDMatrixXSecRW->Write();
  fFDVsNDMatrixTauXSecRW->Write();
                        
  fTrueEnergyNuFlux_ND->Write();
  fTrueEnergyNuFluxRW_ND->Write();
  fTrueEnergyNuFlux_FD->Write();
  fTrueEnergyNuFluxRW_FD->Write();
  fTrueEnergyCCFlux_ND->Write();
  fTrueEnergyCCFluxRW_ND->Write();
  fTrueEnergyCCFlux_FD->Write();
  fTrueEnergyCCFluxRW_FD->Write();
  savefile->Close();
  if (savefile) {delete savefile; savefile = 0;}
}

//____________________________________________________________________72
void NuFluxHelper::FillNonMMHelpers(const NuFluxChain* fluxChain,
                                    const Int_t flukaEntry) const
{
  //Get the fluka variables I'm going to need
  NuParticleType_t Ntype = fluxChain->NuType(flukaEntry);
  Int_t tptype = fluxChain->ParentParticleType(flukaEntry);
  Double_t tpx = fluxChain->TargetParentPX(flukaEntry);
  Double_t tpy = fluxChain->TargetParentPY(flukaEntry);
  Double_t tpz = fluxChain->TargetParentPZ(flukaEntry);
  Double_t Nenergyn = fluxChain->NeutrinoNDEnergy(flukaEntry);
  Double_t Nenergyf = fluxChain->NeutrinoFDEnergy(flukaEntry);
  Double_t Nimpwt = fluxChain->NuImportanceWeight(flukaEntry);
  Double_t Nwtnear = fluxChain->NuNDWeight(flukaEntry);
  Double_t Nwtfar = fluxChain->NuFDWeight(flukaEntry);
  
  //Higher-level fluka variables
  Double_t pt = 0;
  Double_t pt2 = tpx*tpx + tpy*tpy;
  if (0 < pt2){pt = sqrt(pt2);}

  //Make NuEvents to get the beam weights into.  
  NuEvent nearEvent;
  nearEvent.reweightVersion = this->ReweightVersion();
  nearEvent.simFlag = SimFlag::kMC;
  nearEvent.beamType = this->BeamType();
  nearEvent.runPeriod = this->RunPeriod();
  nearEvent.neuEnMC = Nenergyn;
  nearEvent.inu=fluxChain->NuParticleAsPDGCode(Ntype);
  nearEvent.tptype = tptype;
  nearEvent.tpz = tpz;
  nearEvent.tpy = tpy;
  nearEvent.tpx = tpx;
  nearEvent.detector = Detector::kNear;
  nearEvent.beamWeightRunI = 1.0;
  nearEvent.beamWeightRunII = 1.0;
  nearEvent.beamWeight = 1.0;
  nearEvent.fluxErr = 1.0;
  nearEvent.fluxErrTotalErrorAfterTuneRunI = 1.0;
  nearEvent.fluxErrTotalErrorAfterTuneRunII = 1.0;
  nearEvent.fluxErrTotalErrorAfterTune = 1.0;

  NuEvent farEvent;
  farEvent.reweightVersion = this->ReweightVersion();
  farEvent.simFlag = SimFlag::kMC;
  farEvent.beamType = this->BeamType();
  farEvent.runPeriod = this->RunPeriod();
  farEvent.neuEnMC = Nenergyf;
  farEvent.inu=fluxChain->NuParticleAsPDGCode(Ntype);
  farEvent.tptype = tptype;
  farEvent.tpz = tpz;
  farEvent.tpy = tpy;
  farEvent.tpx = tpx;
  farEvent.detector = Detector::kFar;
  farEvent.beamWeightRunI = 1.0;
  farEvent.beamWeightRunII = 1.0;
  farEvent.beamWeight = 1.0;
  farEvent.fluxErr = 1.0;
  farEvent.fluxErrTotalErrorAfterTuneRunI = 1.0;
  farEvent.fluxErrTotalErrorAfterTuneRunII = 1.0;
  farEvent.fluxErrTotalErrorAfterTune = 1.0;

  //Get the SKZP weights
  Double_t flux_weight_near = 1.0;
  Double_t flux_weight_far = 1.0;

  if (fdoSKZP){
    //Reweight the NuEvents
    fzarko->GetBeamWeightsOnly(nearEvent);
    fzarko->GetBeamWeightsOnly(farEvent);
//    Not needed now that run period is specified in the event
//    if (SKZPWeightCalculator::kRunI == frunPeriod){
//      {
//      static Int_t messageCounter = 0;
//      if (messageCounter < 5){
//        MSG("NuFluxHelper",Msg::kInfo)
//          << "Reweighting for run 1, beam type " 
//          << BeamType::AsString(this->BeamType()) << endl;
//        ++messageCounter;
//      }
//      }
//      flux_weight_near = nearEvent.beamWeightRunI;
//      flux_weight_far = farEvent.beamWeightRunI;
//      nearEvent.fluxErr = nearEvent.fluxErrTotalErrorAfterTuneRunI;
//      farEvent.fluxErr = farEvent.fluxErrTotalErrorAfterTuneRunI;
//    }
//    else if (SKZPWeightCalculator::kRunII == frunPeriod){
//      {
//      static Int_t messageCounter = 0;
//      if (messageCounter < 5){
//        MSG("NuFluxHelper",Msg::kInfo)
//          << "Reweighting for run 2, beam type " 
//          << BeamType::AsString(this->BeamType()) << endl;
//        ++messageCounter;
//      }
//      }
//      flux_weight_near = nearEvent.beamWeightRunII;
//      flux_weight_far = farEvent.beamWeightRunII;
//      nearEvent.fluxErr = nearEvent.fluxErrTotalErrorAfterTuneRunII;
//      farEvent.fluxErr = farEvent.fluxErrTotalErrorAfterTuneRunII;
//    }
//    else{
//      {
//      static Int_t messageCounter = 0;
//      if (messageCounter < 5){
//        MSG("NuFluxHelper",Msg::kInfo)
//          << "Reweighting for PoT averaged runs 1 & 2, beam type " 
//          << BeamType::AsString(this->BeamType()) << endl;
//        ++messageCounter;
//      }
//      }
//      flux_weight_near = nearEvent.beamWeight;
//      flux_weight_far = farEvent.beamWeight;
//      nearEvent.fluxErr = nearEvent.fluxErrTotalErrorAfterTune;
//      farEvent.fluxErr = farEvent.fluxErrTotalErrorAfterTune;
//    }
  }
  else{
    static Int_t messageCounter = 0;
    if (messageCounter < 5){
      MSG("NuFluxHelper",Msg::kInfo)
        << "No SKZP reweighting" << endl;
      ++messageCounter;
    }
  }
  if (fdoBeamShift){
    {
      static Int_t messageCounter = 0;
      if (messageCounter < 5){
        MSG("NuFluxHelper",Msg::kInfo)
          << "Performing beam shift "
          << fbeamShiftAsSigma << " sigma"
          << endl;
        ++messageCounter;
      }
    }
    flux_weight_near *= 1.0 + fbeamShiftAsSigma*(nearEvent.fluxErr - 1.0);
    flux_weight_far *= 1.0 + fbeamShiftAsSigma*(farEvent.fluxErr - 1.0);
  }
  if (fdoScrapingShift){
    Double_t ppvx = fluxChain->ParentProdVtxX(flukaEntry);
    Double_t ppvy = fluxChain->ParentProdVtxY(flukaEntry);
    Double_t ppvz = fluxChain->ParentProdVtxZ(flukaEntry);
    Float_t ppvr = sqrt(ppvx*ppvx + ppvy*ppvy);
    Float_t Znom = 52.06; // -187.06 for HE
    if (BeamType::kL250z200i == this->BeamType()){
      Znom = -187.06;
    }
    else if (BeamType::kL010z185i == this->BeamType()){
      Znom = 52.06;
    }
    else{
      Znom = 52.06;
    }
    
    Bool_t shouldIShift = true; 
    if (ppvr < 1.65 && TMath::Abs(ppvz - Znom) < 0.1){
      shouldIShift = false;
    }
    if (TMath::Abs(ppvr - 1.51) < 0.11 && ppvz < Znom){
      shouldIShift = false; // Target Side
    }
    // Chase (z < 4500) or Decay Pipe (z > 4500)
    if (shouldIShift){
      flux_weight_near *= fscrapingScaleFactor;
      flux_weight_far *= fscrapingScaleFactor;
    }
  }
  
  //What's the cross-section of my neutrino?
  //Only needed for filling XSec-reweighted plots
  Double_t xsec_nd = this->CrossSection(Ntype,Nenergyn);
  Double_t xsec_fd = this->CrossSection(Ntype,Nenergyf);
  if(xsec_nd<=0.0) xsec_nd = 0.0;
  if(xsec_fd<=0.0) xsec_fd = 0.0;
  
  //Fill the beam matrix without MM reweighting or cross-sections:
  //(Does this reduce down to the N/F method?)
  fFDVsNDMatrix->Fill(Nenergyn,
                      Nenergyf,
                      Nimpwt*Nwtfar*flux_weight_far);
  
  //Fill the beam matrix without MM reweighting but with cross-sections:
  fFDVsNDMatrixXSec->Fill(Nenergyn,
                          Nenergyf,
                          Nimpwt*Nwtfar*xsec_fd*flux_weight_far);
  
  //Needed to normalise fFDVsNDMatrix:
  fTrueEnergyNuFlux_ND->Fill(Nenergyn,
                             Nimpwt*Nwtnear*flux_weight_near);
  
  //Not needed:
  fTrueEnergyNuFlux_FD->Fill(Nenergyf,
                             Nimpwt*Nwtfar*flux_weight_far);   
  
  //Needed to normalise fFDVsNDMatrixXSec:
  fTrueEnergyCCFlux_ND->Fill(Nenergyn,
                             Nimpwt*Nwtnear*xsec_nd*flux_weight_near);
  
  //Not needed:
  fTrueEnergyCCFlux_FD->Fill(Nenergyf,
                             Nimpwt*Nwtfar*xsec_fd*flux_weight_far);

  return;
}

//____________________________________________________________________72
const Double_t
NuFluxHelper::CrossSection(const NuParticleType_t particle,
                           const Double_t energy) const
{
  if (NuParticle::kNuMu == particle){
    return (Double_t) fxSecGraphNuMuCC->Eval(energy,0,"");
  }
  if (NuParticle::kNuMuBar == particle){
    return (Double_t) fxSecGraphNuMuBarCC->Eval(energy,0,"");
  }
  return 0.0;
}

//____________________________________________________________________72
const Double_t
NuFluxHelper::TauCrossSection(const NuParticleType_t particle,
                              const Double_t energy) const
{
  //The cross section if the mus change to taus en route
  if (NuParticle::kNuMu == particle){
    return (Double_t) fxSecGraphNuTauCC->Eval(energy,0,"");
  }
  if (NuParticle::kNuMuBar == particle){
    return (Double_t) fxSecGraphNuTauBarCC->Eval(energy,0,"");
  }
  return 0.0;
}

//____________________________________________________________________72
const TVector3 NuFluxHelper::
ConvertNDToBeamCoOrdinates(const TVector3& ndCoordinates) const
{
  //Give it a point in ND co-ordinates.
  //Returns the same point in beam co-ordinates.
  
  // beam center (0,0) in detector coord system:
  // (148.844,13.97)
  // z offset is 104000
  // beamdydz = -0.0581
  Double_t beamCenterXDet = 148.28*Munits::cm;
  Double_t beamCenterYDet = 23.84*Munits::cm;
  Double_t beamAngle = 3.323155*TMath::DegToRad();
  Double_t detZ0InBeamCoOrds = 103648.837*Munits::cm;
  
  Double_t beamX = ndCoordinates.X() - beamCenterXDet;
  Double_t beamY =
    (ndCoordinates.Y()-beamCenterYDet)
    *TMath::Cos(beamAngle) +
    ndCoordinates.Z()*TMath::Sin(beamAngle);
  Double_t beamZ =
    detZ0InBeamCoOrds +
    (ndCoordinates.Z()*TMath::Cos(beamAngle) - 
     ndCoordinates.Y()*TMath::Sin(beamAngle));
  //N.B. beamz should have a term beamCenterYDet*TMath::Sin(3.23155)
  //assuming the 104000 cm is measured along the beam axis.
  //But this term's small.

  TVector3 beamCoordinates(beamX,
                           beamY,
                           beamZ);
  return beamCoordinates;
}

//____________________________________________________________________72
const Bool_t NuFluxHelper::IsNDFiducial(const TVector3& ndPoint) const
{
  //create a NuCuts object
  static NuCuts nuCuts;
  
  //give fake values for u,v,plane,r,releaseType
  //since they are not needed for kCC0325Std
  return nuCuts.IsInFidVol(ndPoint.X(),ndPoint.Y(),ndPoint.Z(),
                           0,0,
                           0,0,
                           Detector::kNear,NuCuts::kCC0325Std,
                           0,SimFlag::kMC);
  
  //Bool_t NuCuts::IsInFidVol(Float_t x,Float_t y,Float_t z,
  //              Float_t u,Float_t v,
  //              Int_t planeTrkVtx,Float_t rTrkVtx,
  //              Int_t detector,Int_t anaVersion,
  //              Int_t releaseType,Int_t simFlag)
  
  //return nuCuts.IsInFidVolNDCC0325Std(ndPoint.X(),
  //                        ndPoint.Y(),
  //                        ndPoint.Z());
}

//____________________________________________________________________72
void NuFluxHelper::FillMMHelpers(const NuFluxChain* fluxChain,
                                 const Int_t flukaEntry) const
{
  //Get the fluka variables I'm going to need
  NuParticleType_t Ntype = fluxChain->NuType(flukaEntry);
  Int_t tptype = fluxChain->ParentParticleType(flukaEntry);
  Double_t tpx = fluxChain->TargetParentPX(flukaEntry);
  Double_t tpy = fluxChain->TargetParentPY(flukaEntry);
  Double_t tpz = fluxChain->TargetParentPZ(flukaEntry);
  Double_t Nimpwt = fluxChain->NuImportanceWeight(flukaEntry);
  
  //Higher-level fluka variables
  Double_t pt = 0;
  Double_t pt2 = tpx*tpx + tpy*tpy;
  if (0 < pt2){pt = sqrt(pt2);}
  
  //Loop over the neutrinos 10 times to spread out importance weights
  for(int subLoop=0;subLoop<10;subLoop++){
    
    //Make a vector ready to hold random detector points
    TVector3 vDetPoint(0,0,0);
    
    //Keep creating random points until one falls within the
    //fiducial volume.
    do {
      vDetPoint.SetX(frandom3->Uniform(-300,300)*Munits::cm);
      vDetPoint.SetY(frandom3->Uniform(-300,300)*Munits::cm);
      vDetPoint.SetZ(frandom3->Uniform(0,600)*Munits::cm);
    } while (!this->IsNDFiducial(vDetPoint));
    
    TVector3 vBeamPoint = this->ConvertNDToBeamCoOrdinates(vDetPoint);
    
    Double_t X = vBeamPoint.X()/Munits::cm;
    Double_t Y = vBeamPoint.Y()/Munits::cm;
    Double_t Z = vBeamPoint.Z()/Munits::cm;
    
    //Get the ND weight and energy:
    ArrAy newvals_nd = this->NuWte(X,Y,Z,fluxChain,flukaEntry);
    //get new fd weight and energy:
    ArrAy newvals_fd = this->NuWte(0,0,73534000.,fluxChain,flukaEntry);
    
    //Calculate SKZP parameters for the new neutrino energies
    //Make NuEvents to get the beam weights into.  
    NuEvent nearEvent;
    nearEvent.reweightVersion = this->ReweightVersion();
    nearEvent.simFlag = SimFlag::kMC;
    nearEvent.beamType = this->BeamType();
    nearEvent.runPeriod = this->RunPeriod();
    nearEvent.neuEnMC = newvals_nd.new_ene;
    nearEvent.inu=fluxChain->NuParticleAsPDGCode(Ntype);
    nearEvent.tptype = tptype;
    nearEvent.tpz = tpz;
    nearEvent.tpy = tpy;
    nearEvent.tpx = tpx;
    nearEvent.detector = Detector::kNear;
    nearEvent.beamWeightRunI = 1.0;
    nearEvent.beamWeightRunII = 1.0;
    nearEvent.beamWeight = 1.0;
    nearEvent.fluxErr = 1.0;
    nearEvent.fluxErrTotalErrorAfterTuneRunI = 1.0;
    nearEvent.fluxErrTotalErrorAfterTuneRunII = 1.0;
    nearEvent.fluxErrTotalErrorAfterTune = 1.0;
    
    NuEvent farEvent;
    farEvent.reweightVersion = this->ReweightVersion();
    farEvent.simFlag = SimFlag::kMC;
    farEvent.beamType = this->BeamType();
    farEvent.runPeriod = this->RunPeriod();
    farEvent.neuEnMC = newvals_fd.new_ene;
    farEvent.inu=fluxChain->NuParticleAsPDGCode(Ntype);
    farEvent.iaction=1;//CC
    farEvent.tptype = tptype;
    farEvent.tpz = tpz;
    farEvent.tpy = tpy;
    farEvent.tpx = tpx;
    farEvent.detector = Detector::kFar;
    farEvent.beamWeightRunI = 1.0;
    farEvent.beamWeightRunII = 1.0;
    farEvent.beamWeight = 1.0;
    farEvent.fluxErr = 1.0;
    farEvent.fluxErrTotalErrorAfterTuneRunI = 1.0;
    farEvent.fluxErrTotalErrorAfterTuneRunII = 1.0;
    farEvent.fluxErrTotalErrorAfterTune = 1.0;
    
    //Get the SKZP weights
    Double_t flux_weight_near = 1.0;
    Double_t flux_weight_far = 1.0;
    
    if (fdoSKZP){
      //Reweight the NuEvents
      fzarko->GetBeamWeightsOnly(nearEvent);
      fzarko->GetBeamWeightsOnly(farEvent);
      
      // Debugging code -- you shouldn't see this!
      if (farEvent.beamWeightRunI != 1.0) {
        cerr << "Beam weights are wrong -- beware!!!" << endl;
        assert(false);
      }
      
      flux_weight_near = nearEvent.beamWeight;
      flux_weight_far = farEvent.beamWeight;
      nearEvent.fluxErr = nearEvent.fluxErrTotalErrorAfterTune;
      farEvent.fluxErr = farEvent.fluxErrTotalErrorAfterTune;
    }
    if (fdoBeamShift){
      {
        static Int_t messageCounter = 0;
        if (messageCounter < 5){
          MSG("NuFluxHelper",Msg::kInfo)
            << "Performing beam shift "
            << fbeamShiftAsSigma << " sigma"
            << endl;
          ++messageCounter;
        }
      }
      flux_weight_near *= 1.0 + fbeamShiftAsSigma*(nearEvent.fluxErr - 1.0);
      flux_weight_far *= 1.0 + fbeamShiftAsSigma*(farEvent.fluxErr - 1.0);
    }
    if (fdoScrapingShift){
      Double_t ppvx = fluxChain->ParentProdVtxX(flukaEntry);
      Double_t ppvy = fluxChain->ParentProdVtxY(flukaEntry);
      Double_t ppvz = fluxChain->ParentProdVtxZ(flukaEntry);
      Float_t ppvr = sqrt(ppvx*ppvx + ppvy*ppvy);
      Float_t Znom = 52.06; // -187.06 for HE
      if (BeamType::kL250z200i == this->BeamType()){
        Znom = -187.06;
      }
      else if (BeamType::kL010z185i == this->BeamType()){
        Znom = 52.06;
      }
      else{
        Znom = 52.06;
      }
      
      Bool_t shouldIShift = true; 
      if (ppvr < 1.65 && TMath::Abs(ppvz - Znom) < 0.1){
        shouldIShift = false;
      }
      if (TMath::Abs(ppvr - 1.51) < 0.11 && ppvz < Znom){
        shouldIShift = false; // Target Side
      }
      // Chase (z < 4500) or Decay Pipe (z > 4500)
      if (shouldIShift){
        flux_weight_near *= fscrapingScaleFactor;
        flux_weight_far *= fscrapingScaleFactor;
      }
    }
    
    //Calculate cross-sections for the new neutrino energies
    Double_t xsec_nd = this->CrossSection(Ntype,newvals_nd.new_ene);
    if(xsec_nd<=0.0) xsec_nd = 0.0; 
    Double_t xsec_fd = this->CrossSection(Ntype,newvals_fd.new_ene);
    if(xsec_fd<=0.0) xsec_fd = 0.0; 
    //Tau cross section:
    Double_t tauXsec_fd = this->TauCrossSection(Ntype,newvals_fd.new_ene);
    if (tauXsec_fd<=0.0) tauXsec_fd = 0.0;
    
    //Fill the beam matrix:
    fFDVsNDMatrixXSecRW->Fill(newvals_nd.new_ene,
                              newvals_fd.new_ene,
                              Nimpwt*newvals_fd.new_weight*
                              xsec_fd*flux_weight_far);

    //Fill the numu->nutau beam matrix:
    fFDVsNDMatrixTauXSecRW->Fill(newvals_nd.new_ene,
                                 newvals_fd.new_ene,
                                 Nimpwt*newvals_fd.new_weight*
                                 tauXsec_fd*flux_weight_far);
    
    
    //Fill the beam matrix without cross-sections:
    fFDVsNDMatrixRW->Fill(newvals_nd.new_ene,
                          newvals_fd.new_ene,
                          Nimpwt*newvals_fd.new_weight*
                          flux_weight_far);
    
    //Needed to normalise fFDVsNDMatrixRW:
    fTrueEnergyNuFluxRW_ND->Fill(newvals_nd.new_ene,
                                 Nimpwt*newvals_nd.new_weight*
                                 flux_weight_near);
    //Not needed:
    fTrueEnergyNuFluxRW_FD->Fill(newvals_fd.new_ene,
                                 Nimpwt*newvals_fd.new_weight*
                                 flux_weight_far);
    
    //Needed to normalise fFDVsNDMatrixXSecRW and fFDVsNDMatrixTauXSecRW:
    fTrueEnergyCCFluxRW_ND->Fill(newvals_nd.new_ene,
                                 Nimpwt*newvals_nd.new_weight*
                                 xsec_nd*flux_weight_near);
    //Not needed:
    fTrueEnergyCCFluxRW_FD->Fill(newvals_fd.new_ene,
                                 Nimpwt*newvals_fd.new_weight*
                                 xsec_fd*flux_weight_far);
       
  } //End x10 loop over fluka neutrinos
}

//____________________________________________________________________72
void NuFluxHelper::NormaliseHistos() const
{
  MSG("NuFluxHelper",Msg::kInfo)
    <<"Normalising histograms" << endl;

  //Loop over the bins of the beam matrices:
  for(int i=1;i<=fNNDTrueEBins;i++){
    for(int j=1;j<=fNFDTrueEBins;j++){
      //Renormalise the un-MM-reweighted, un-XSec-reweighted beam matrix
      if(fTrueEnergyNuFlux_ND->GetBinContent(i)>0 && 
         fFDVsNDMatrix->GetBinContent(i,j)>0 ) {
        Float_t error = TMath::Power(fFDVsNDMatrix->GetBinError(i,j) / 
                                     fFDVsNDMatrix->GetBinContent(i,j),2);
        error = TMath::Sqrt(error);
        fFDVsNDMatrix->SetBinContent(i,j,
                                     fFDVsNDMatrix->GetBinContent(i,j)/
                                     fTrueEnergyNuFlux_ND->GetBinContent(i));
        fFDVsNDMatrix->SetBinError(i,j,error * 
                                   fFDVsNDMatrix->GetBinContent(i,j));
      }
      else {
        fFDVsNDMatrix->SetBinContent(i,j,0.);
        fFDVsNDMatrix->SetBinError(i,j,0.);
      }
      
      //Renormalise the MM-reweighted, un-XSec-reweighted beam matrix
      if( fTrueEnergyNuFluxRW_ND->GetBinContent(i)>0 && 
          fFDVsNDMatrixRW->GetBinContent(i,j)>0 ) {
        Float_t error = TMath::Power(fFDVsNDMatrixRW->GetBinError(i,j) / 
                                     fFDVsNDMatrixRW->GetBinContent(i,j),2);
        error = TMath::Sqrt(error);
        fFDVsNDMatrixRW->SetBinContent(i,j,
                                       fFDVsNDMatrixRW->GetBinContent(i,j)/
                                       fTrueEnergyNuFluxRW_ND->GetBinContent(i));
        fFDVsNDMatrixRW->SetBinError(i,j,error * 
                                     fFDVsNDMatrixRW->GetBinContent(i,j));
      }
      else {
        fFDVsNDMatrixRW->SetBinContent(i,j,0.);
        fFDVsNDMatrixRW->SetBinError(i,j,0.);
      }

      //Renormalise the un-MM-reweighted, XSec-reweighted beam matrix
      if( fTrueEnergyCCFlux_ND->GetBinContent(i)>0 &&
          fFDVsNDMatrixXSec->GetBinContent(i,j)>0 ) {
        Float_t error = TMath::Power(fFDVsNDMatrixXSec->GetBinError(i,j) / 
                                     fFDVsNDMatrixXSec->GetBinContent(i,j),2);
        error = TMath::Sqrt(error);
        fFDVsNDMatrixXSec->SetBinContent(i,j,
                                         fFDVsNDMatrixXSec->GetBinContent(i,j)/
                                         fTrueEnergyCCFlux_ND->GetBinContent(i));
        fFDVsNDMatrixXSec->SetBinError(i,j,error * 
                                       fFDVsNDMatrixXSec->GetBinContent(i,j));
      }
      else {
        fFDVsNDMatrixXSec->SetBinContent(i,j,0.);
        fFDVsNDMatrixXSec->SetBinError(i,j,0.);
      }

      //Renormalise the MM-reweighted, XSec-reweighted beam matrix
      if( fTrueEnergyCCFluxRW_ND->GetBinContent(i)>0 &&
          fFDVsNDMatrixXSecRW->GetBinContent(i,j)>0 ) {
        Float_t error = TMath::Power(fFDVsNDMatrixXSecRW->GetBinError(i,j) / 
                                     fFDVsNDMatrixXSecRW->GetBinContent(i,j),2);
        error = TMath::Sqrt(error);
        fFDVsNDMatrixXSecRW->SetBinContent(i,j,
                                           fFDVsNDMatrixXSecRW->GetBinContent(i,j)/
                                           fTrueEnergyCCFluxRW_ND->GetBinContent(i));
        fFDVsNDMatrixXSecRW->SetBinError(i,j,error *
                                         fFDVsNDMatrixXSecRW->GetBinContent(i,j));
      }
      else {
        fFDVsNDMatrixXSecRW->SetBinContent(i,j,0.);
        fFDVsNDMatrixXSecRW->SetBinError(i,j,0.);
      }

      //Renormalise the MM-reweighted, tauXSec-reweighted beam matrix
      if( fTrueEnergyCCFluxRW_ND->GetBinContent(i)>0 &&
          fFDVsNDMatrixTauXSecRW->GetBinContent(i,j)>0 ) {
        Float_t error = TMath::Power(fFDVsNDMatrixTauXSecRW->GetBinError(i,j) / 
                                     fFDVsNDMatrixTauXSecRW->GetBinContent(i,j),2);
        error = TMath::Sqrt(error);
        fFDVsNDMatrixTauXSecRW->SetBinContent(i,j,
                                              fFDVsNDMatrixTauXSecRW->GetBinContent(i,j)/
                                              fTrueEnergyCCFluxRW_ND->GetBinContent(i));
        fFDVsNDMatrixTauXSecRW->SetBinError(i,j,error *
                                            fFDVsNDMatrixTauXSecRW->GetBinContent(i,j));
      }
      else {
        fFDVsNDMatrixTauXSecRW->SetBinContent(i,j,0.);
        fFDVsNDMatrixTauXSecRW->SetBinError(i,j,0.);
      }

    } //End loop over matrix columns
  } //End loop over matrix rows
  return;
}

//____________________________________________________________________72
const ArrAy NuFluxHelper::NuWte(const Double_t XDET,
                                const Double_t YDET,
                                const Double_t ZDET,
                                const NuFluxChain* fluxChain,
                                const Int_t flukaEntry) const
{
  ArrAy   newvars,empty;
  empty.new_ene   =-1;
  empty.new_weight=-1.;
  
  Int_t Ntype = fluxChain->NuType(flukaEntry);
  Double_t Vx = fluxChain->ParentDecayVtxX(flukaEntry)/Munits::cm;
  Float_t Vy = fluxChain->ParentDecayVtxY(flukaEntry)/Munits::cm;
  Float_t Vz = fluxChain->ParentDecayVtxZ(flukaEntry)/Munits::cm;
  Float_t pdpx = fluxChain->ParentPX(flukaEntry);
  Float_t pdpy = fluxChain->ParentPY(flukaEntry);
  Float_t pdpz = fluxChain->ParentPZ(flukaEntry);
  Float_t ppdxdz = fluxChain->ParentProdDXDZ(flukaEntry);
  Float_t ppdydz = fluxChain->ParentProdDYDZ(flukaEntry);
  Float_t pppz = fluxChain->ParentProdPZ(flukaEntry);
  Float_t ppenergy = fluxChain->ParentProdEnergy(flukaEntry);
  Int_t ptype = fluxChain->ParentType(flukaEntry);
  Float_t muparpx = fluxChain->MuonParentPX(flukaEntry);
  Float_t muparpy = fluxChain->MuonParentPY(flukaEntry);
  Float_t muparpz = fluxChain->MuonParentPZ(flukaEntry);
  Float_t mupare = fluxChain->MuonParentEnergy(flukaEntry);
  Float_t Necm = fluxChain->NuECofM(flukaEntry);
  
  //    -----------------fixed parameters
  Double_t pimass, kmass, k0mass, mumass, omegamass;
  pimass=0.13957; kmass=0.49368;
  k0mass=0.49767; mumass=0.105658389;
  omegamass=1.67245;
  Int_t nue, nuebar, numu, numubar, muplus, muminus;
  nue=53; nuebar=52; numu=56; numubar=55;
  muplus=5;  muminus=6;
  //      set to flux per 1 meter radius 
  Double_t RDET =100.;
  //   
  
  
//    -----------------ADAMO style variables
      Double_t Neutrino_type;
      Double_t Vertex_x,           Vertex_y,          Vertex_z;
      Double_t Neutrino_parentpx,  Neutrino_parentpy, Neutrino_parentpz;
      Double_t Particle_dxdz,      Particle_dydz,     Particle_pz;
      Double_t Particle_energy,    Particle_type;
      Double_t Muparent_px,        Muparent_py,       Muparent_pz;
      Double_t Muparent_energy;
      Double_t Neutrino_ecm;

//    -----------------rest of local variables
      Double_t eneu, WGT;
      double ENUZR, gamma, beta[3], SANGDET, RAD, EMRAT;
      double beta_mag, gamma_sqr;
      double parent_energy, parent_mass, parentp, partial;
      double costh_pardet, THETA_pardet,costh;
      Double_t P_dcm_nu[4], P_nu[3], P_pcm_mp[3];
      Double_t P_pcm, wt_ratio, xnu;
//======================================================================

//    switch to ADAMO variables to maintain compatibility with internal GNUMI
      Neutrino_type = Ntype;
      Vertex_x = Vx;
      Vertex_y = Vy;
      Vertex_z = Vz;
      Neutrino_parentpx = pdpx;
      Neutrino_parentpy = pdpy;
      Neutrino_parentpz = pdpz;
      Particle_dxdz   = ppdxdz;
      Particle_dydz   = ppdydz;
      Particle_pz     = pppz;
      Particle_energy = ppenergy;
      Particle_type   = ptype;
      Muparent_px = muparpx;
      Muparent_py = muparpy;
      Muparent_pz = muparpz;
      Muparent_energy = mupare;
      Neutrino_ecm = Necm;

//    >=t gamma of muon parent particle at point of decay to neutrino

      if(Particle_type==8||Particle_type==9)        parent_mass = pimass;
      else if(Particle_type==11||Particle_type==12) parent_mass = kmass;
      else if(Particle_type==10||Particle_type==16) parent_mass = k0mass;
      else if(Particle_type==5||Particle_type==6)   parent_mass = mumass;
      else if(Particle_type==24||Particle_type==32) parent_mass = omegamass;
      else{
        cout << "NUWEIGHT unknown particle type STOP" << endl;
        return empty;
      }
      parent_energy = sqrt(double(Neutrino_parentpx)*double(Neutrino_parentpx)
       + double(Neutrino_parentpy)*double(Neutrino_parentpy)
       + double(Neutrino_parentpz)*double(Neutrino_parentpz) + parent_mass*parent_mass);
      gamma = parent_energy / parent_mass;
      gamma_sqr = gamma*gamma;
      beta_mag  = sqrt( (gamma_sqr-1.)/gamma_sqr );

//    >=t the neutrino energy in the parent decay cm

//    :::ARBITRARY LORENTZ TRANSFORM
//    :::     EtP = PCX * BGX + PCY * BGY + PCZ * BGZ
//    :::     E  = GA * E//+ EtP
//    :::     PtE = EtP / (GA + 1e0) + EC
//    :::     PX = PCX + BGX * PtE
//    :::     PY = PCY + BGY * PtE
//    :::     PZ = PCZ + BGZ * PtE
//    :::     P  = SQRT (PX * PX + PY * PY + PZ * PZ)
      ENUZR = Neutrino_ecm;

//    >=t angle from parent line of flight to detector in lab frame

      RAD =pow((double(XDET)-double(Vertex_x)),2)+pow((double(YDET)-double(Vertex_y)),2)
        +pow((double(ZDET)-double(Vertex_z)),2);
      RAD = sqrt(RAD);
      parentp =  double(Neutrino_parentpx)*double(Neutrino_parentpx) +
       double(Neutrino_parentpy)*double(Neutrino_parentpy) + double(Neutrino_parentpz)* double(Neutrino_parentpz);
      parentp = sqrt(parentp);
      costh_pardet = ( Neutrino_parentpx*(double(XDET)-Vertex_x)
       + Neutrino_parentpy*(double(YDET)-Vertex_y)
       + Neutrino_parentpz*(double(ZDET)-Vertex_z) )
       / ( parentp * RAD );
      if(fabs(costh_pardet)>1.){
        if(costh_pardet<0) costh_pardet = -1.; 
        if(costh_pardet>0) costh_pardet =  1;
      }
      THETA_pardet = acos(costh_pardet);

//    weighted neutrino energy in lab, approx, good for small theta
//I'VE WORKED THIS OUT AND I DON'T THINK IT'S APPROXIMATE:
//I THINK IT'S EXACT --JE
      EMRAT = 1. / (gamma * (1. - beta_mag * cos(THETA_pardet)));
      eneu=EMRAT*ENUZR;

//    >=t solid angle/4pi for detector element
      SANGDET = (   RDET*RDET/((ZDET-Vertex_z)*(ZDET-Vertex_z)))/(4.0   );

//    weight for solid angle and lorentz boost
      WGT = SANGDET * (EMRAT*EMRAT);
//    weight weighted by approx neutrino cross section
//    WGTE=0.67*ENEU*WGT
      
//    done for all except polarized muon decay, ==================
//    in which case need to modify weight
//    (Warning: This section may need more double precision)
      if(Particle_type==muplus||Particle_type==muminus) {

//      boost new neutrino to mu decay cm
//      boost new neutrino to mu decay cm
        beta[0] = Neutrino_parentpx / parent_energy;
        beta[1] = Neutrino_parentpy / parent_energy;
        beta[2] = Neutrino_parentpz / parent_energy;
        P_nu[0] = (XDET-Vertex_x)*eneu/RAD;
        P_nu[1] = (YDET-Vertex_y)*eneu/RAD;
        P_nu[2] = (ZDET-Vertex_z)*eneu/RAD;
        partial =gamma*(beta[0]*P_nu[0]+beta[1]*P_nu[1]+beta[2]*P_nu[2]);
        partial = eneu - partial /(gamma+1.);
        P_dcm_nu[0] = P_nu[0] - beta[0]*gamma*partial;
        P_dcm_nu[1] = P_nu[1] - beta[1]*gamma*partial;
        P_dcm_nu[2] = P_nu[2] - beta[2]*gamma*partial;
        P_dcm_nu[3] = sqrt(P_dcm_nu[0]*P_dcm_nu[0]+P_dcm_nu[1]*+P_dcm_nu[1]+P_dcm_nu[2]*P_dcm_nu[2]);

//      boost parent of mu to mu production cm
        gamma = Particle_energy/parent_mass;
        beta[0] = Particle_dxdz*Particle_pz/Particle_energy;
        beta[1] = Particle_dydz*Particle_pz/Particle_energy;
        beta[2] =               Particle_pz/Particle_energy;
        partial = gamma * ( beta[0]*Muparent_px + 
         beta[1]*Muparent_py + beta[2]*Muparent_pz);
        partial = Muparent_energy - partial /(gamma+1.);
        P_pcm_mp[0] = Muparent_px - beta[0]*gamma*partial;
        P_pcm_mp[1] = Muparent_py - beta[1]*gamma*partial;
        P_pcm_mp[2] = Muparent_pz - beta[2]*gamma*partial;
        P_pcm = sqrt(P_pcm_mp[0]*P_pcm_mp[0] + P_pcm_mp[1]* P_pcm_mp[1] + P_pcm_mp[2]* P_pcm_mp[2]);

//      cal//new  decay angle w.r.t. (anti)spin direction
        if(P_dcm_nu[3]*P_pcm!=0) {
          costh  = ( P_dcm_nu[0]*P_pcm_mp[0] + P_dcm_nu[1]*P_pcm_mp[1]
                     + P_dcm_nu[2]*P_pcm_mp[2] ) / (P_dcm_nu[3]*P_pcm);
        }
        else costh = 0;
        if(fabs(costh)>=1.) {
         if(costh<0.) costh = -1. ; 
         if(costh>0.) costh=  1.  ;
        } 
//      cal//relative weight due to angle difference
        if( Neutrino_type==nue || Neutrino_type==nuebar ) {
          wt_ratio = 1.-costh;
        }  
        else if( Neutrino_type==numu|| Neutrino_type==numubar) {
          xnu = 2. * ENUZR / mumass;
          wt_ratio = ( (3.-2.*xnu) - (1.-2.*xnu)*costh ) / (3.-2.*xnu);
         }
        else {
          cout << "NUWEIGHT bad neutrino type" << endl;
          return empty;
        }

        WGT = WGT * wt_ratio;
      }
      Double_t ene_wei[2];
      ene_wei[0]=eneu;
      ene_wei[1]=WGT;
      newvars.new_ene    = eneu;
      newvars.new_weight = WGT;
      return newvars;
}

---