NCEventInfo.cxx

Go to the documentation of this file.

#include "NCEventInfo.h"

#include <cassert>

#include "AnalysisNtuples/ANtpAnalysisInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpDefaultValue.h"
#include "AnalysisNtuples/ANtpEventInfo.h"
#include "AnalysisNtuples/ANtpEventInfoNC.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpShowerInfo.h"
#include "AnalysisNtuples/ANtpShowerInfoNC.h"
#include "AnalysisNtuples/ANtpTrackInfoAtm.h"
#include "AnalysisNtuples/ANtpTrackInfoNC.h"
#include "AnalysisNtuples/ANtpTruthInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoAtm.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"

#include "Conventions/Detector.h"
#include "Conventions/Munits.h"

#include "DataUtil/EnergyCorrections.h"

#include "MCReweight/MCReweight.h"
#include "MCReweight/SKZPWeightCalculator.h"

#include "MessageService/MsgService.h"

#include "Validity/VldTimeStamp.h"
#include "Validity/VldContext.h"


#include "TChain.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH2F.h"
#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"

#include "NCUtils/NCOscProb.h"
#include "NCUtils/NCType.h"
#include "NCUtils/NCRunUtil.h"
#include "NCUtils/NCUtility.h"
using NC::Utility::SQR;

using namespace EnergyCorrections;

CVSID("$Id: NCEventInfo.cxx,v 1.49 2010/01/19 17:02:22 rodriges Exp $");

static const int kQE = 0;
static const int kRES = 1;
static const int kDIS = 2;

// define static member
NCEventInfo::INukeHists NCEventInfo::sfINukeHists;

NCEventInfo::INukeHists::INukeHists()
{
  //save the current working directory
  TDirectory* tmp = gDirectory;

  //  fCuts = new NCAnalysisCutsNC();
  //  fZbeam = new Zbeam();
  //  fZfluk = new Zfluk();

  //  fSKZPWeight = new SKZPWeightCalculator(beamWeightConfig,true);
  //  fZbeam->SetReweightConfig(beamWeightConfig);

  // ---RPL 05/09/06: instantiate the mcReweight object now, and give
  // it a (neugen) weight calculator
  // Back this off; Brian says it should be done in the script.
  //MCReweight &mcReweight = MCReweight::Instance();
  //fNeugenWeightCal = new NeugenWeightCalculator();
  //mcReweight.AddWeightCalculator(fNeugenWeightCal);
  // ---

  //gotta load shower energy weighting histograms from Mike K.
  TString scaleFromName = "integrated_smearing_histos_default.root";
  TString scaleToNameMinus = "integrated_smearing_histos_1501.root";
  TString scaleToNamePlus = "integrated_smearing_histos_1500.root";
  TString offsetName = "integrated_smearing_histos_offset.root";

  TString topDir="MCReweight/data";
  TString base=getenv("SRT_PRIVATE_CONTEXT");
  if(base!="" && base!="."){
    // check if directory exists in SRT_PRIVATE_CONTEXT
    TString path = base + "/" + topDir;
    void *dir_ptr = gSystem->OpenDirectory(path);
    if(!dir_ptr) base=getenv("SRT_PUBLIC_CONTEXT");
  }
  // if it doesn't exist then use SRT_PUBLIC_CONTEXT
  else base=getenv("SRT_PUBLIC_CONTEXT");

  if(base==""){
    MSG("NCUtils",Msg::kFatal) << "No SRT_PUBLIC_CONTEXT set" << endl;
    assert(false);
  }
  topDir = base+ "/" + topDir + "/";

  MSG("NCUtils",Msg::kDebug)
    << "Zbeam reading files from: " << topDir << endl;
  TFile scaleFileFrom(topDir+scaleFromName);
  TFile scaleFileToMinus(topDir+scaleToNameMinus);
  TFile scaleFileToPlus(topDir+scaleToNamePlus);
  TFile offsetFile(topDir+offsetName);

  TString names[] = {"h_qe", "h_res", "h_dis"};

  TH2F *histFrom = dynamic_cast<TH2F *>(scaleFileFrom.Get(names[kQE]));
  TH2F *histToMinus = dynamic_cast<TH2F *>(scaleFileToMinus.Get(names[kQE]));
  TH2F *histToPlus = dynamic_cast<TH2F *>(scaleFileToPlus.Get(names[kQE]));
  TH2F *hist = dynamic_cast<TH2F *>(offsetFile.Get(names[kQE]));

  for(int i = kQE; i < kDIS+1; ++i){

    histFrom = dynamic_cast<TH2F *>(scaleFileFrom.Get(names[i]));
    histToMinus = dynamic_cast<TH2F *>(scaleFileToMinus.Get(names[i]));
    histToPlus = dynamic_cast<TH2F *>(scaleFileToPlus.Get(names[i]));
    hist = dynamic_cast<TH2F *>(offsetFile.Get(names[i]));

    MSG("NCEventInfo", Msg::kDebug)
      << histFrom->GetName() << " "
      << histToMinus->GetName() << " "
      << histToPlus->GetName() << " "
      << hist->GetName() << endl;

    from.push_back(histFrom);
    toMinus.push_back(histToMinus);
    toPlus.push_back(histToPlus);
    offset.push_back(hist);

    from[i]->SetDirectory(0);
    toMinus[i]->SetDirectory(0);
    toPlus[i]->SetDirectory(0);
    offset[i]->SetDirectory(0);

    MSG("NCEventInfo", Msg::kDebug)
      << from[i]->GetName() << " "
      << toMinus[i]->GetName() << " "
      << toPlus[i]->GetName() << " "
      << offset[i]->GetName() << endl;

  }//end loop to fill histograms

  gDirectory = tmp;
}

//---------------------------------------------------------------------
NCEventInfo::INukeHists::~INukeHists()
{
  for(unsigned int n = 0; n < offset.size(); ++n) delete offset[n];
  for(unsigned int n = 0; n < from.size(); ++n) delete from[n];
  for(unsigned int n = 0; n < toMinus.size(); ++n) delete toMinus[n];
  for(unsigned int n = 0; n < toPlus.size(); ++n) delete toPlus[n];
}

//---------------------------------------------------------------------
NCEventInfo::NCEventInfo(const char* beamWeightConfig)
  : analysis(0), beam(0), event(0), header(0),
    reco_nom(0), reco(0), shower(0), track(0), truth(0),
    fSKZPWeight(0),
    fBeamWeightConfig(beamWeightConfig)
{
}

//---------------------------------------------------------------------
SKZPWeightCalculator* NCEventInfo::GetSKZPCalc()
{
  if(!fSKZPWeight)
    fSKZPWeight = new SKZPWeightCalculator(fBeamWeightConfig, true);

  const ReleaseType::Release_t rel
    = ReleaseType::StringToType(header->softVersion);

  if(fBeamWeightConfig == "PiMinus_CedarDaikon"){
    if(!ReleaseType::IsCedar(rel)){
      MAXMSG("NCEventInfo", Msg::kError, 10) << endl
        << "*******************************************************" << endl
        << "* APPLYING CEDARDAIKON SKZP TO NON CEDAR FILE         *" << endl
        << "* FIX BEFORE YOU DO A REAL ANALYSIS.                  *" << endl
        << "*******************************************************" << endl;
    }
    if(ReleaseType::IsMC(rel) && !ReleaseType::IsDogwood(rel)){
      MAXMSG("NCEventInfo", Msg::kError, 10) << endl
        << "*******************************************************" << endl
        << "* APPLYING CEDARDAIKON SKZP TO NON DOGWOOD FILE       *" << endl
        << "* FIX BEFORE YOU DO A REAL ANALYSIS.                  *" << endl
        << "*******************************************************" << endl;
    }
  }

  if(fBeamWeightConfig == "Dogwood1_Daikon07" && !ReleaseType::IsDogwood(rel)){
    MAXMSG("NCEventInfo", Msg::kError, 10) << endl
      << "*******************************************************" << endl
      << "* APPLYING DOGWOODDAIKON07 SKZP TO NON DOGWOOD FILE   *" << endl
      << "* FIX BEFORE YOU DO A REAL ANALYSIS.                  *" << endl
      << "*******************************************************" << endl;
  }
  

  return fSKZPWeight;
}

//---------------------------------------------------------------------
void NCEventInfo::DeepCopy(const NCEventInfo* evt)
{
#define COPYFIELD(fld, typ) \
  if(evt->fld){             \
    if(!fld) fld = new typ; \
    *fld = *evt->fld;       \
  }                         \
  else{                     \
    if(fld) delete fld;     \
    fld = 0;                \
  }

  COPYFIELD(analysis, ANtpAnalysisInfo);
  COPYFIELD(beam,     ANtpBeamInfo);
  COPYFIELD(event,    ANtpEventInfoNC);
  COPYFIELD(header,   ANtpHeaderInfo);
  COPYFIELD(reco,     ANtpRecoInfo);
  COPYFIELD(reco_nom, ANtpRecoInfo);
  COPYFIELD(shower,   ANtpShowerInfoNC);
  COPYFIELD(track,    ANtpTrackInfoNC);
  COPYFIELD(truth,    ANtpTruthInfoBeam);

#undef COPYFIELD
}

//---------------------------------------------------------------------
double NCEventInfo::GetEventVertex(double& x, double& y, double& z) const
{
  // choice of whether to use Event or Track vertex should follow what is in
  //NCAnalysisCuts::InBeamFiducialVolumeOx(), otherwise the vertex that is being
  //cut and the vertex that is being stored are not the same

  assert(header);
  assert(event && track && shower);

  // If we have reco_nom then it's analysis time and we shouldn't
  // be redoing the calculation but merely using the variables we
  // calculated before.
  assert(!reco_nom);

  if(event->tracks > 0 &&
     track->planes - shower->planes > 0){
    x = track->vtxX;
    y = track->vtxY;
    z = track->vtxZ  - NCType::kTrackVtxAdjustment;
  }
  else{
    x = event->vtxX;
    y = event->vtxY;
    z = event->vtxZ;
  }

  // Correct the coordinates into beam coordinates for the purposes of
  // calculating this return value only.
  if(header->detector == int(Detector::kNear)){
    const double x1 = x - NCType::kNDBeamCenterX;
    const double y1 = y - NCType::kNDBeamCenterY - NCType::kNDBeamAngle*z;
    return TMath::Sqrt(x1*x1 + y1*y1);
  }

  return TMath::Sqrt(x*x + y*y);
}


//---------------------------------------------------------------------
double NCEventInfo::GetEventEnergy(CandShowerHandle::ShowerType_t showerType,
                                   bool stoppingBeamMuon) const
{
  double showerEnergy = 0.;
  double trackEnergy = 0.;

  if(event->showers > 0)
    showerEnergy = GetShowerEnergy(showerType);
  if(event->tracks > 0)
    trackEnergy = GetTrackEnergy(stoppingBeamMuon);//fCuts->IsStoppingBeamMuon());


  if(showerType == CandShowerHandle::kWtNC) return showerEnergy;

  return trackEnergy + showerEnergy;
}


//---------------------------------------------------------------------
double NCEventInfo::GetTrackEnergy(bool contained) const
{
  assert(header);

  const SimFlag::SimFlag_t sim = SimFlag::SimFlag_t(header->dataType);

  using namespace ReleaseType;
  const Release_t rt = StringToType(header->softVersion);

  const Detector::Detector_t det = Detector::Detector_t(header->detector);

  //get a validity context for the energy corrections
  VldContext vldc(det, sim, GetTimestamp());

  bool fromCurve = false;
  bool fromRange = false;

  double inputMomentum;

  if(contained || track->fitMomentum < -9000) {
    if(track->rangeMomentum > 0){
      fromRange = true;
      inputMomentum = track->rangeMomentum;
    }
    else{
      fromCurve = true;
      inputMomentum = track->fitMomentum;
    }
  }
  else{
    fromCurve = true;
    inputMomentum = track->fitMomentum;
  }

  assert(fromCurve || fromRange);
  assert(!(fromCurve && fromRange));

  // required to shut gcc up in optimized mode
  // the asserts above guarantee that one of the branches below will be taken
  double trackMomentum = 0;

  if(fromCurve)
    trackMomentum = FullyCorrectSignedMomentumFromCurvature(inputMomentum,
                                                            vldc, rt);
  if(fromRange)
    trackMomentum = FullyCorrectMomentumFromRange(inputMomentum, vldc, rt);

  return TMath::Sqrt(SQR(trackMomentum) + SQR(NCType::kMuMassGeV));
}

//---------------------------------------------------------------------
double NCEventInfo::doubleintpow(double a, int b) const
{
  double ret = 1;
  while(b--) ret *= a;
  return ret;
}

//---------------------------------------------------------------------

double NCEventInfo::MasakiStyleCorrectionImp(double logE, const double pars[6]) const
{
  // logE is log_10(shower energy)
  // pars are the Chebyshev coefficients
  double ret = 2.-(pars[0]+
                   pars[1]*logE+
                   pars[2]*(2.*doubleintpow(logE, 2)-1.)+
                   pars[3]*(4.*doubleintpow(logE, 3)-3.*logE)+
                   pars[4]*(8.*doubleintpow(logE, 4)-8.*doubleintpow(logE, 2)+1.)+
                   pars[5]*(16.*doubleintpow(logE, 5)-
                            20.*doubleintpow(logE, 3)+5.*logE));
  ret = min(ret, 1.3);
  ret = max(ret, 0.7);
  return ret;
}

//---------------------------------------------------------------------
double NCEventInfo::MasakiStyleCorrectionCedarPhyLinfix(double E) const
{
  // Masaki-style energy corrections appropriate for the
  // cedar_phy_linfix sample. See docdb 5753

  const ReleaseType::Release_t rel
    = ReleaseType::StringToType(header->softVersion);
  if(!ReleaseType::IsCedar(rel)){
    MAXMSG("NCEventInfo", Msg::kError, 10) << endl
      << "***********************************************************" << endl
      << "* APPLYING CEDAR-PHY ENERGY CORRECTIONS TO NON CEDAR FILE *" << endl
      << "* FIX BEFORE YOU DO A REAL ANALYSIS.                      *" << endl
      << "***********************************************************" << endl;
  }
  if(ReleaseType::IsMC(rel) && !ReleaseType::IsDogwood(rel)){
    MAXMSG("NCEventInfo", Msg::kError, 10) << endl
      << "***********************************************************" << endl
      << "* APPLYING ENERGY CORRECTIONS TO NON DOGWOOD FILE         *" << endl
      << "* FIX BEFORE YOU DO A REAL ANALYSIS.                      *" << endl
      << "***********************************************************" << endl;
  }

  if(E==0) return 1.;

  E=max(E, 0.2);

  // Different corrections for near/far, low/high energy (E<>10 GeV)

  const double parsNDLo[6]={ 8.78922e-01,  2.01774e-01, -1.91411e-01,
                             1.45199e-01, -1.02706e-01 , 3.67004e-02 };
  const double parsNDHi[6]={ 1.02466,      2.76558e-01, -6.60202e-01,
                             4.36941e-01, -1.21538e-01,  0.0120561   };
  const double parsFDLo[6]={ 8.49032e-01,  2.81519e-01, -2.48210e-01,
                             1.73536e-01, -9.64045e-02,  3.11845e-02 };
  const double parsFDHi[6]={ 9.66207e-01,  3.04379e-01, -5.45190e-01,
                             3.49389e-01, -9.28317e-02,  0.00870293  };

  const double* parsToUse;

  if(header->detector == (int)Detector::kNear)
    parsToUse= (E<=10) ? parsNDLo : parsNDHi;
  else if(header->detector == (int)Detector::kFar)
    parsToUse= (E<=10) ? parsFDLo : parsFDHi;
  else assert (0 && "DISASTER!!!!!! Unknown detector");

  return MasakiStyleCorrectionImp(log10(E), parsToUse);
}

//---------------------------------------------------------------------
double NCEventInfo::
GetShowerEnergy(CandShowerHandle::ShowerType_t showerType) const
{
  assert(header);
  assert(shower);

  // If we have reco_nom then it's analysis time and we shouldn't
  // be redoing the calculation but merely using the variables we
  // calculated before.
  assert(!reco_nom);

  const SimFlag::SimFlag_t sim = SimFlag::SimFlag_t(header->dataType);

  using namespace ReleaseType;
  const Release_t rt = StringToType(header->softVersion);

  const Detector::Detector_t det = Detector::Detector_t(header->detector);

  //get a validity context for the energy corrections
  VldContext vldc(det, sim, GetTimestamp());

  double showerEnergy = ANtpDefaultValue::kDouble;

  switch(showerType){
  case CandShowerHandle::kCC:
    showerEnergy = shower->linearCCGeV;
    break;
  case CandShowerHandle::kWtCC:
    showerEnergy = shower->deweightCCGeV;
    break;
  case CandShowerHandle::kNC:
    showerEnergy = shower->linearNCGeV;
    break;
  case CandShowerHandle::kWtNC:
    showerEnergy = shower->deweightNCGeV;
    break;
  case CandShowerHandle::kEM:
    assert(0 && "Not implemented");
  }

  if(showerEnergy > 0){
    showerEnergy = FullyCorrectShowerEnergy(showerEnergy,
                                            showerType, vldc, rt);

    // Do post-hoc MEU corrections for cedar_phy_linfix MC. Taken from
    // docdb 5340 for the ND, and an email I have from Alex (forwarded
    // from Greg) for the FD
    //
    // These constants are only for cedar_phy_linfix, but there's no way
    // to tell that from the release type so cedar is the best we can test for
    if(sim == SimFlag::kMC && ReleaseType::IsCedar(rt)){
      if(det == Detector::kNear)     showerEnergy *= 1.007667;
      else if(det == Detector::kFar) showerEnergy *= 1.0018;
      else MSG("NCEventInfo", Msg::kWarning) << "Don't know detector: "
                                             << det << endl;
    }//MC

    // Only apply the Masaki-style corrections to CC showers
    if(showerType != CandShowerHandle::kCC) return showerEnergy;
    return (showerEnergy*MasakiStyleCorrectionCedarPhyLinfix(showerEnergy));
  }

  return showerEnergy;
}


//---------------------------------------------------------------------------
//will find MODBYRS3 weight for MC version carrot or earlier even when all
//values of useParameters are set to false.
double NCEventInfo::FindNeugenWeight(const bool* useParameters,
                                     const std::vector<double>& adjustedValues,
                                     bool useMODBYRS3, int overrideMODBYRS)
{
  assert(truth);

  TString neugen = "neugen:";

  MCReweight &mcReweight = MCReweight::Instance();

  int nucleus = 1; //default

  //figure out the nucleus
  if(int(truth->atomicNumber == 1))  nucleus = 0; // free
  else if(int(truth->atomicNumber == 6))  nucleus = 274; // oxygen
  else if(int(truth->atomicNumber == 8))  nucleus = 284; // carbon
  else if(int(truth->atomicNumber == 26)) nucleus = 372; // iron

  if(nucleus == 1) return 1;

  MCEventInfo ei;
  ei.nuE           = truth->nuEnergy;
  ei.nuPx          = truth->nuDCosX*truth->nuEnergy;
  ei.nuPy          = truth->nuDCosY*truth->nuEnergy;
  ei.nuPz          = truth->nuDCosZ*truth->nuEnergy;
  ei.tarE          = truth->targetEnergy;
  ei.tarPx         = truth->targetPX;
  ei.tarPy         = truth->targetPY;
  ei.tarPz         = truth->targetPZ;
  ei.y             = truth->hadronicY;
  ei.x             = truth->bjorkenX;
  ei.q2            = truth->q2;
  ei.w2            = truth->w2;
  ei.iaction       = truth->interactionType;
  ei.inu           = truth->nuFlavor;
  ei.iresonance    = truth->resonanceCode;
  ei.initial_state = truth->initialState;
  ei.nucleus       = nucleus;
  ei.had_fs        = truth->hadronicFinalState;

  NuParent* np = 0;

  Registry fReweightConfigRegistry;
  fReweightConfigRegistry.UnLockValues();
  Registry defaultConfigRegistry;

  const char* mcString = (header->softVersion).Data();
  ReleaseType::Release_t mcType=ReleaseType::StringToType(mcString);


  //Only apply MODBYRS3 for Carrot or earlier
  //Note: According to discussion with Hugh G. and
  //Tricia V., Carrot was actually generated with MODBYRS2
  //and then found to agree better with data using MODBYRS3.
  //Anyone using this code on carrot will get unweighted carrot if
  //useMODBYRS3=false and reweighted carrot if flag useMODBYRS3=true

  //overriding default version values for MODBYRS
  if(overrideMODBYRS != 0)
  {
     fReweightConfigRegistry.Set("neugen:config_name", "MODBYRS");
     fReweightConfigRegistry.Set("neugen:config_no", overrideMODBYRS);
  }
  //no override
  else
  {
     if(useMODBYRS3 && ReleaseType::IsCarrot(mcType) ){
        fReweightConfigRegistry.Set("neugen:config_name", "MODBYRS");
        fReweightConfigRegistry.Set("neugen:config_no", 3);
     }
     else if(ReleaseType::IsDaikon(mcType)){
        fReweightConfigRegistry.Set("neugen:config_name", "MODBYRS");
        fReweightConfigRegistry.Set("neugen:config_no", 4);
     }
  }
  /*
  if((int) useParameters.size() < NCType::kNumNeugenParameters){
    MSG("NCEventInfo", Msg::kWarning) << "use parameters not same size"
                                          << " as neugen parameters - return"
                                          << " weight of 1." << endl;
    return 1.;
  }
  */

  //set the configuration for the reweight based on the percent change
  double change = 1.;
  for(int i = 0; i < (int)adjustedValues.size(); ++i){
    change = adjustedValues[i];

    //for now the reweighting doesnt work right - see comment below
    //make sure that the defaults for all unchanged parameters are the
    //current defaults
    if(!useParameters[i]
       && i < NCType::kDISFACT
       && !(i == NCType::kma_qe && useParameters[NCType::kCCMA])
       && !(i == NCType::kma_res && useParameters[NCType::kCCMA])
       && !(i == NCType::kkno_r112 && useParameters[NCType::kkno_r112122])
       && !(i == NCType::kkno_r122 && useParameters[NCType::kkno_r112122])
       && !(i == NCType::kkno_r113 && useParameters[NCType::kkno_r113123])
       && !(i == NCType::kkno_r123 && useParameters[NCType::kkno_r113123])
       && !(i == NCType::kkno_r212 && useParameters[NCType::kkno_r212222])
       && !(i == NCType::kkno_r222 && useParameters[NCType::kkno_r212222])
       && !(i == NCType::kkno_r213 && useParameters[NCType::kkno_r213223])
       && !(i == NCType::kkno_r223 && useParameters[NCType::kkno_r213223])) {

       if(ReleaseType::IsDaikon(mcType))
       {
          fReweightConfigRegistry.Set(neugen+NCType::kParams[i].name,
                                      NCType::kParams[i].defaultVal);
       }
    }

    if(useParameters[i] && i < NCType::kDISFACT){
      fReweightConfigRegistry.Set(neugen+NCType::kParams[i].name, change);
    }
    else if(useParameters[i] && i == NCType::kDISFACT){
      //here adjusted values is the fractional change for the parameters,
      //multiply it by the default values
      fReweightConfigRegistry.Set("neugen:kno_r112",
                                  change*NCType::kParams[NCType::kkno_r112].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r122",
                                  change*NCType::kParams[NCType::kkno_r122].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r132",
                                  change*NCType::kParams[NCType::kkno_r132].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r142",
                                  change*NCType::kParams[NCType::kkno_r142].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r113",
                                  change*NCType::kParams[NCType::kkno_r113].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r123",
                                  change*NCType::kParams[NCType::kkno_r123].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r133",
                                  change*NCType::kParams[NCType::kkno_r133].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r143",
                                  change*NCType::kParams[NCType::kkno_r143].defaultVal);
    }
    else if(useParameters[i] && i == NCType::kCCMA){
      fReweightConfigRegistry.Set("neugen:ma_qe",
                                  change*NCType::kParams[NCType::kma_qe].defaultVal);
      fReweightConfigRegistry.Set("neugen:ma_res",
                                  change*NCType::kParams[NCType::kma_res].defaultVal);
    }
    else if(useParameters[i] && i == NCType::kkno_r112122){
      fReweightConfigRegistry.Set("neugen:kno_r112",
                                  change*NCType::kParams[NCType::kkno_r112].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r122",
                                  change*NCType::kParams[NCType::kkno_r122].defaultVal);
    }
    else if(useParameters[i] && i == NCType::kkno_r113123){
      fReweightConfigRegistry.Set("neugen:kno_r113",
                                  change*NCType::kParams[NCType::kkno_r113].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r123",
                                  change*NCType::kParams[NCType::kkno_r123].defaultVal);
    }
    else if(useParameters[i] && i == NCType::kkno_r212222){
      fReweightConfigRegistry.Set("neugen:kno_r212",
                                  change*NCType::kParams[NCType::kkno_r212].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r222",
                                  change*NCType::kParams[NCType::kkno_r222].defaultVal);
    }
    else if(useParameters[i] && i == NCType::kkno_r213223){
      fReweightConfigRegistry.Set("neugen:kno_r213",
                                  change*NCType::kParams[NCType::kkno_r213].defaultVal);
      fReweightConfigRegistry.Set("neugen:kno_r223",
                                  change*NCType::kParams[NCType::kkno_r223].defaultVal);
    }

  }//end loop over adjusted values

  mcReweight.ResetAllReweightConfigs();

  double weight = mcReweight.ComputeWeight(&ei,np,&fReweightConfigRegistry);

  //reweighting is in strange state right now - it returns the weight
  //with respect to parameters as of MODBYRS2 - Daikon is MODBYRS4
  //so have to get weight from MODBYRS2 defaults to MODBYRS4 defaults first
  double defaultWeight = 1.;
  if(ReleaseType::IsDaikon(mcType)){
    defaultConfigRegistry.Set("neugen:config_name", "MODBYRS");
    defaultConfigRegistry.Set("neugen:config_no", 4);

    //use defaults for all parameters up to the combinations
    for(int i = 0; i < NCType::kDISFACT; ++i)
      defaultConfigRegistry.Set(neugen+NCType::kParams[i].name,
                                NCType::kParams[i].defaultVal);

    mcReweight.ResetAllReweightConfigs();

    defaultWeight = mcReweight.ComputeWeight(&ei,np,&defaultConfigRegistry);

  }//end if daikon

  //divide reweight by default weight
  //cout << "oldw = " << weight << endl;
  if(TMath::Abs(defaultWeight) > 0.) weight /= defaultWeight;
  //cout << "neww = " << weight << endl;
  //cout << "neww*w = " << weight * reco->weight << endl;
  if( TMath::IsNaN(weight) ){
    MSG("NCEventInfo", Msg::kWarning)
      << "mcreweight IsNaN "
      <<  truth->nuEnergy
      << " " << truth->interactionType
      << " " << truth->resonanceCode
      << " " << truth->nuFlavor
      << " " << truth->initialState
      << " " << truth->hadronicFinalState
      << " " << nucleus
      << endl;
    weight = 1.;
  }

  if(weight != 1.){
    MAXMSG("NCEventInfo", Msg::kDebug, 100)
      << "mcreweight to " << weight << endl
      << "values into neugen" << endl
      << "event:nu_en "
      << truth->nuEnergy << endl
      << "event:nu_px "
      << truth->nuDCosX*truth->nuEnergy << endl
      << "event:nu_py "
      << truth->nuDCosY*truth->nuEnergy << endl
      << "event:nu_pz "
      << truth->nuDCosZ*truth->nuEnergy << endl
      << "event:tar_en "
      << truth->targetEnergy << endl
      << "event:tar_px "
      << truth->targetPX << endl
      << "event:tar_py "
      << truth->targetPY << endl
      << "event:tar_pz "
      << truth->targetPZ << endl
      << "event:x "
      << truth->bjorkenX << endl
      << "event:y "
      << truth->hadronicY << endl
      << "event:q2 " << truth->q2 << endl
      << "event:w2 " << truth->w2 << endl
      << "event:iaction "
      << truth->interactionType << endl
      << "event:inu "
      << truth->nuFlavor << endl
      << "event:iresonance "
      << truth->resonanceCode << endl
      << "event:initial_state "
      << truth->initialState << endl
      << "event:nucleus " << nucleus << endl
      << "event:had_fs "
      << truth->hadronicFinalState
      << endl << endl;
  }

  return weight;
}


//---------------------------------------------------------------------------
double NCEventInfo::FindCrossSectionWeight(const bool* useParameters,
                                           const std::vector<double>& adjustedValues) const
{
  double weight = 1.;
  double trueE = truth->nuEnergy;

  if(useParameters[NCType::kNCCrossSection]
     && truth->interactionType == NCType::kNC){
    if(trueE <= 10.)
      weight = 1. + 0.2*adjustedValues[NCType::kNCCrossSection];
    else if(trueE < 100.)
      weight = 1. + 0.2*(1. - 0.0111*(trueE-10.))*adjustedValues[NCType::kNCCrossSection];
    else
      weight = 1.;
  }

  if(useParameters[NCType::kNuBarCrossSection]
     && truth->interactionType == NCType::kCC
     && truth->nuFlavor < 0)
    weight *= 1. + adjustedValues[NCType::kNuBarCrossSection];

  return weight;
}


//----------------------------------------------------------------------
//calls FindNeugenWeight without any parameters set to true.  FindNeugenWeight
//returns the MODBYRS3 weight in that case if the release is carrot or
//earlier, otherwise returns 1.
double NCEventInfo::FindMODBYRSWeight()
{
  vector<double> adjustedValues;
  bool useParameters[NCType::kNumNeugenParameters] = {false,};

  for(int i = 0; i < NCType::kNumNeugenParameters; ++i){
    adjustedValues.push_back(NCType::kParams[i].defaultVal);
  }

  return FindNeugenWeight(useParameters, adjustedValues);
}


//---------------------------------------------------------------------------
//this version of the method uses the shiny new functionality
//in MCReweight to get the #'s from the db
double NCEventInfo::FindMEGAFitWeight(int beamType, NC::RunUtil::ERunType runType)
{
  //check if there is a beam object loaded, if so use it to get the beam type,
  //otherwise use the passed in int.
  int beamt = BeamType::ToZarko((BeamType::BeamType_t)(beamType)); //Use SKZP beam encoding
  if(beam) beamt = BeamType::ToZarko((BeamType::BeamType_t)(beam->beamType)); //Use SKZP beam encoding

  SKZPWeightCalculator::RunPeriod_t rp;

  // Bah, converting between conventions
  switch(runType){
  case NC::RunUtil::kRunI:
    rp=SKZPWeightCalculator::kRunI;
    break;
  case NC::RunUtil::kRunII:
    rp=SKZPWeightCalculator::kRunII;
    break;
  case NC::RunUtil::kRunIII:
    rp=SKZPWeightCalculator::kRunIII;
    break;
  default:
    assert(0 && "Unknown run period");
    break;
  }

  double pt = TMath::Sqrt(SQR(truth->targetParentPX)
                          +SQR(truth->targetParentPY));
  double emuNew = 0.;
  double eshwNew = 0.;
  double beamweight = 0.;
  double detweight = 0.;
  double weightbeam = GetSKZPCalc()->GetWeight(header->detector,
                                               beamt,
                                               truth->targetParentType,
                                               pt,
                                               1.*truth->targetParentPZ,
                                               truth->interactionType,
                                               truth->nuEnergy,
                                               truth->nuFlavor,
                                               0.,
                                               0.,
                                               emuNew,
                                               eshwNew,
                                               beamweight,
                                               detweight,
                                               rp);

  //only do the MODBYRS weighting if using carrot
  double weightMODBYRS = 1.;
  const char* mcString = (header->softVersion).Data();
  ReleaseType::Release_t mcType = ReleaseType::StringToType(mcString);
  if( ReleaseType::IsCarrot(mcType) ) weightMODBYRS = FindMODBYRSWeight();

  MAXMSG("NCEventInfo", Msg::kDebug, 10)
    << "MODBYRS3 = " << weightMODBYRS
    << " beam = " << weightbeam
    << endl;

  return weightMODBYRS*weightbeam;
}


//---------------------------------------------------------------------------
double NCEventInfo::FindMEGAFitWeightUncertainty(int beamType,
                                                 NC::RunUtil::ERunType runType,
                                                 double sigma,
                                                 SKZPWeightCalculator* skzpcalc)
{
  assert(header && truth);

  SKZPWeightCalculator::RunPeriod_t rp = SKZPWeightCalculator::kRunI;
  if(runType == NC::RunUtil::kRunII) rp = SKZPWeightCalculator::kRunII;

  return skzpcalc->GetFluxError(header->detector,
                                BeamType::ToZarko((BeamType::BeamType_t)(beamType)), //use SKZP enum
                                truth->nuFlavor,
                                truth->nuEnergy,
                                SKZPWeightCalculator::kTotalErrorAfterTune,
                                sigma);
}


//---------------------------------------------------------------------------
double NCEventInfo::
CalculateAbsoluteHadronicShiftScale(double trueShowerEnergy,
                                    double sigma) const
{
  const double calibrationError = 5.7;

  assert(trueShowerEnergy >= 0);

  double modellingError = -1;
  if(trueShowerEnergy < 0.5) modellingError = 8.2;
  else if(trueShowerEnergy < 10) modellingError = 2.7+3.7*TMath::Exp(-.25*trueShowerEnergy);
  else modellingError = 3;
  assert(modellingError > 0);

  const double totalError = TMath::Sqrt(SQR(calibrationError)+
                                        SQR(modellingError));
  return 1+.01*totalError*sigma;
}


//---------------------------------------------------------------------------

void NCEventInfo::ShiftEnergies(const bool* useParameters,
                                const vector<double>& adjustedValues)
{
  assert(analysis);
  assert(reco_nom);
  //  assert(useParameters.size() == adjustedValues.size());

  // Reset the energies in the reco object to their versions in
  // reco_nom to ensure that this function is idempotent
  reco->muEnergy=reco_nom->muEnergy;
  reco->showerEnergyNC=reco_nom->showerEnergyNC;
  reco->showerEnergyCC=reco_nom->showerEnergyCC;
  reco->nuEnergyNC=reco_nom->nuEnergyNC;
  reco->nuEnergyCC=reco_nom->nuEnergyCC;

  const double trackEnergy_nom=reco_nom->muEnergy;

  if(TMath::Abs(trackEnergy_nom) > NCType::kMuMassGeV){
    double trackMomentum = TMath::Sqrt(SQR(trackEnergy_nom) - SQR(NCType::kMuMassGeV));

    if(useParameters[NCType::kTrackEnergy])
      trackMomentum *= 1. + adjustedValues[NCType::kTrackEnergy];

    reco->muEnergy = TMath::Sqrt(SQR(trackMomentum) + SQR(NCType::kMuMassGeV));
  }

  double showerEnergyNC_nom = reco_nom->showerEnergyNC;
  double showerEnergyCC_nom = reco_nom->showerEnergyCC;

  // Only try to shift the shower if there is one
  if(! (ANtpDefaultValue::IsDefault(showerEnergyNC_nom) &&
        ANtpDefaultValue::IsDefault(showerEnergyCC_nom)) ){
    //old way of doing it by scaling shower energy
    if(useParameters[NCType::kShowerEnergy] &&
       !useParameters[NCType::kAbsoluteHadronicCalibration]) {
      const double scale= 1. + adjustedValues[NCType::kShowerEnergy];

      reco->showerEnergyNC*=scale;
      reco->showerEnergyCC*=scale;
    }

    if(useParameters[NCType::kAbsoluteHadronicCalibration]){
      const double trueShowerEnergy = truth->showerEnergy;
      const double sigma = adjustedValues[NCType::kAbsoluteHadronicCalibration];
      const double scale = CalculateAbsoluteHadronicShiftScale(trueShowerEnergy, sigma);
      reco->showerEnergyNC *= scale;
      reco->showerEnergyCC *= scale;
    }

    if(useParameters[NCType::kRelativeHadronicCalibration] &&
       header->detector == Detector::kFar){
      const double scale= 1. + adjustedValues[NCType::kRelativeHadronicCalibration];
      
      reco->showerEnergyNC *= scale;
      reco->showerEnergyCC *= scale;
    }
    assert(reco->showerEnergyNC >= 0);
    assert(reco->showerEnergyCC >= 0);
  }

  //new way using mike's histograms

  int offset = kQE;
  int from = kQE;
  int to = kQE;
  SelectINukeHist(offset, from, to);

  // if(useParameters[1] && adjustedValues[1] < 0.)
  //   showerEnergy = ReweightInuke(fINukeFrom[from], fINukeToMinus[to]);
  // else if(useParameters[1] && adjustedValues[1] > 0.)
  //   showerEnergy = ReweightInuke(fINukeFrom[from], fINukeToPlus[to]);

  //Shower offset must be greater than 0 to use.
  if(useParameters[NCType::kShowerEnergyOffset]
     && TMath::Abs(adjustedValues[NCType::kShowerEnergyOffset])>0.0)
    // This is wrong, but will make it compile. If you get here,
    // SimulateShowerOffset will abort anyway
    reco->showerEnergyNC  = SimulateShowerOffset(sfINukeHists.offset[offset],
                                                 adjustedValues[NCType::kShowerEnergyOffset]);

  reco->nuEnergyNC = reco->showerEnergyNC;
  reco->nuEnergyCC = reco->muEnergy + reco->showerEnergyCC;

  if(analysis->isNC>0){ 
    reco->showerEnergy=reco->showerEnergyNC;
    reco->nuEnergy=reco->nuEnergyNC;
  }

  if(analysis->isCC>0){ 
    reco->showerEnergy=reco->showerEnergyCC;
    reco->nuEnergy=reco->nuEnergyCC;
  }
}

//---------------------------------------------------------------------------

double NCEventInfo::FindRecoWeight(const bool* useParameters,
                                   const vector<double>& adjustedValues) const
{
  double weight = 1;

  const double showerEnergyNC=reco->showerEnergyNC;

  //adjust the weight for the event based on the parameters to fit
  if(useParameters[NCType::kNormalization]
     && header->dataType==SimFlag::kMC
     && header->detector==Detector::kFar)
    weight *= 1.+adjustedValues[NCType::kNormalization];

  if(useParameters[NCType::kNCBackground]
     && truth->interactionType < 1
     && analysis->isCC > 0)
    weight *= 1.+adjustedValues[NCType::kNCBackground];

  //CC background is real CC's, not less than 50% complete CC's.
  //9/4/07: BJR - >50% complete is still chosen as CC after the data cleaning
  if(useParameters[NCType::kCCBackground]
     && truth->interactionType > 0
     && TMath::Abs(truth->nuFlavor) == 14
     && analysis->isNC > 0
     //     && truth->eventCompleteness > 0.5
     ) {
    weight *= 1.+adjustedValues[NCType::kCCBackground];
  }
  if(useParameters[NCType::kPIDCut]
     && analysis->separationParameterPAN < adjustedValues[NCType::kPIDCut])
    weight = 0.;

  //new Values for the Low Completeness. It should not be used alone, so it can be dropped in a decond time
  if(useParameters[NCType::kLowCompleteness]
     && analysis->isNC > 0
     && truth->trueVisibleE &&
     ((showerEnergyNC/truth->trueVisibleE)<0.3)){
    if(showerEnergyNC <= 0.5)
      weight *= 1.+0.13*adjustedValues[NCType::kLowCompleteness];
    else if (showerEnergyNC > 0.5 && showerEnergyNC <=1.0)
      weight *= 1.+0.83*adjustedValues[NCType::kLowCompleteness];
    else weight *= 1.+adjustedValues[NCType::kLowCompleteness];
  }

  //NCFarCleanNoise +-7.8% if Ereco<0.5 GeV +-1.5% if 0.5<EReco<0.75GeV
  if(useParameters[NCType::kNCFarCleanNoise]
     && header->detector==Detector::kFar
     && analysis->isNC > 0){
    if(showerEnergyNC <= 0.5)
      weight  *= 1. + 0.078*adjustedValues[NCType::kNCFarCleanNoise];
    else if (showerEnergyNC > 0.5 && showerEnergyNC <=0.75)
      weight  *= 1. + 0.015*adjustedValues[NCType::kNCFarCleanNoise];
    else weight *= 1.;
  }

  //NCFarCleanCR +-0.0161*e^(-Ereco/6.96)
  if(useParameters[NCType::kNCFarCleanCR]
     && header->detector==Detector::kFar
     && analysis->isNC > 0){
    weight  *= 1. + (0.0161*TMath::Exp(-showerEnergyNC)/6.96)*adjustedValues[NCType::kNCFarCleanCR];
  }

  //Systematics for isSimpleCutsClean
  //OVERALL cleaning, cuts+Poorly reconstructed events
  //+-6%  Ereco< 1GeV
  //+-2%  E <1.5 GeV
  //+-1%  E>1.5 GeV
  if(useParameters[NCType::kNCNearClean] && header->detector==Detector::kNear &&
     analysis->isNC > 0){
    const double cleaningSigmas=adjustedValues[NCType::kNCNearClean];
    if(showerEnergyNC <= 1.){
      //cout << "a. energy=" << showerEnergyNC << " weightFactor=" 
      //     << 1. + 0.06*cleaningSigmas << endl;
      weight *= 1. + 0.06*cleaningSigmas;
    }
    else if(showerEnergyNC > 1. && showerEnergyNC <=1.5){
      //cout << "b. energy=" << showerEnergyNC << " weightFactor=" 
      //     << 1. + 0.02*cleaningSigmas << endl;
      weight *= 1. + 0.02*cleaningSigmas;
    }
    else{
      //cout << "c. energy=" << showerEnergyNC << " weightFactor=" 
      //     << 1. + 0.01*cleaningSigmas << endl;
      weight *= 1. + 0.01*cleaningSigmas;
    }

  }

  //NCCleanRunDiff survey  +-4% if     Ereco<1 GeV in ND
  if(useParameters[NCType::kNCCleanRunDiff]
     && header->detector==Detector::kNear
     && analysis->isNC > 0){
    if(showerEnergyNC <= 1.0)
      weight  *= 1. + 0.04*adjustedValues[NCType::kNCCleanRunDiff];
    else weight *= 1.;
  }


  //NCRunDiff survey  +-4% if     Ereco<=1 GeV
  //                  +-2% if  1 < Ereco < 5 GeV
  if(useParameters[NCType::kNCRunDiff]
       && analysis->isNC > 0){
    if(showerEnergyNC <= 1.0)
      weight  *= 1. + 0.04*adjustedValues[NCType::kNCRunDiff];
    else if(showerEnergyNC > 1.0 && showerEnergyNC < 5.0)
      weight  *= 1. + 0.02*adjustedValues[NCType::kNCRunDiff];
    else weight *= 1.;
  }

  // Messages are slow
  /*
  for(int i = NCType::kTrackEnergy; i < NCType::kNCFarCleanCR; ++i){
    MAXMSG("NCEventInfo", Msg::kDebug, 1) << adjustedValues[i] << " ";
  }
  MAXMSG("NCEventInfo", Msg::kDebug, 1) << endl;

  MAXMSG("NCEventInfo", Msg::kDebug, 30)
    << trackEnergy << "/" << reco->muEnergy << " "
    << showerEnergy << "/" << reco->showerEnergy
    << " " << energy << "/" << reco->nuEnergy
    << " " << weight << endl;
  */

  return weight;
}


//---------------------------------------------------------------------
void NCEventInfo::SelectINukeHist(int& offset, int& from, int& to) const
{
  assert(truth);

  const double trueY = truth->nuEnergy*truth->hadronicY;

  // pick the histograms to use, based on process but with
  // cut on energy. (reason for cut = sparse statistics)
  if(truth->resonanceCode == 1001){
    if(truth->trueVisibleE < 1.0) offset = kQE;
    else if(truth->trueVisibleE < 2.0) offset = kRES;
    else offset = kDIS;

    if(trueY < 1.0){
      from = kQE;
      to = kQE;
    }
    else if(trueY < 2.0){
      from = kRES;
      to = kRES;
    }
    else{
      from = kDIS;
      to = kDIS;
    }
  }//end if QE
  else if(truth->resonanceCode == 1002
          || truth->resonanceCode == 1004){
    if(truth->trueVisibleE < 2.0) offset = kRES;
    else offset = kDIS;

    if(trueY < 2.0){
      from = kRES;
      to = kRES;
    }
    else{
      from = kDIS;
      to = kDIS;
    }
  }
  else{
    from = kDIS;
    to = kDIS;
    offset = kDIS;
  }
}


//---------------------------------------------------------------------
double NCEventInfo::SimulateShowerOffset(TH2F* offsetHist, double offset) const
{
  // (PAR 2009-11-27) I really don't understand this function, and I
  // don't think we use it anyway, so I'm moving the onus for checking
  // it to any potential future users
  MSG("NCEventInfo", Msg::kFatal) 
    << "SimulateShowerOffset is untested."
    << "If you really want to use it, read the code"
    << "and check it does what you want" << endl;
  // Dummy message to make sure it exits
  MSG("NCEventInfo", Msg::kFatal) << "Exiting" << endl;

  assert(truth && reco_nom);

  MAXMSG("NCEventInfo", Msg::kDebug, 20)
    << "simulate offset "
    << offsetHist->GetName()
    << " " << offset
    << " " << truth->trueVisibleE
    << " " << reco_nom->nuEnergy
    << " " << reco_nom->muEnergy
    << " " << reco_nom->showerEnergy
    << endl;

  // don't do anything to IMD
  if(truth->resonanceCode==1005) return reco_nom->showerEnergy;
  // sufficiently far from zero
  if(truth->trueVisibleE > 5.){

    if(reco_nom->showerEnergy+offset < 0.)
      return 0.;

    return reco_nom->showerEnergy+offset;
  }
  // below zero = zero
  if(truth->trueVisibleE+offset < 0.) return 0.;

  // find original bin in reco and true visible energy
  int iox = offsetHist->GetXaxis()->FindBin(truth->trueVisibleE);
  int ioy = offsetHist->GetYaxis()->FindBin(reco_nom->showerEnergy);
  double porig = offsetHist->GetBinContent(iox,ioy);

  // find shifted bin
  int isx = offsetHist->GetXaxis()->FindBin(truth->trueVisibleE
                                            + offset);
  int isy = 0;
  double pnew = 0;
  for(int i = 1; i <= offsetHist->GetNbinsY(); i++){
    pnew = offsetHist->GetBinContent(isx,i);
    if(pnew >= porig){
      isy = i;
      break;
    }
  }

  // if we have moved bins, sample from a uniform distribution
  // in the new bins
  double eup = offsetHist->GetYaxis()->GetBinUpEdge(isy);
  double elow = offsetHist->GetYaxis()->GetBinLowEdge(isy);
  double E = gRandom->Rndm()*(eup-elow) + elow;

  // if we haven't moved bins in reco
  // do not change the energy
  if(isy == ioy) return reco_nom->showerEnergy;

  return E;
}


// UNUSED
/*
//---------------------------------------------------------------------
double NCEventInfo::ReweightInuke(TH2F* intranukeFromHist,
                                  TH2F* intranukeToHist)
{
  double trueShowerE = truth->nuEnergy*truth->hadronicY;

  // don't do anything to IMD
  if(truth->resonanceCode==1005) return reco->showerEnergy;
  // very high energy showers do not suffer much from intranuke
  // also, they are hard to cope with due to poor statistical precision
  if(trueShowerE > 15.0) return reco->showerEnergy;

  // find original bin in reco and true visible energy
  int iox = intranukeFromHist->GetXaxis()->FindBin(trueShowerE);
  int ioy = intranukeFromHist->GetYaxis()->FindBin(reco->showerEnergy);
  double porig = intranukeFromHist->GetBinContent(iox, ioy);

  // find new bin in reco energy
  int isx = intranukeToHist->GetXaxis()->FindBin(trueShowerE);
  int isy = 0;
  double pnew = 0;
  for(int i = 1; i <= intranukeToHist->GetNbinsY(); i++){
    pnew = intranukeToHist->GetBinContent(isx, i);
    if(pnew >= porig){
      isy = i;
      break;
    }
  }

  // if we have moved bins, sample from a uniform distribution
  // in the new bins
  double eup = intranukeToHist->GetYaxis()->GetBinUpEdge(isy);
  double elow = intranukeToHist->GetYaxis()->GetBinLowEdge(isy);
  double E = gRandom->Rndm()*(eup - elow) + elow;

  // if we haven't moved bins in reco
  // do not change the energy
  if(isy == ioy) return reco->showerEnergy;

  return E;
}
*/


//----------------------------------------------------------------------
double NCEventInfo::FindTransitionProbability(const NC::OscProb::OscPars* pars,
                                              bool useNuE) const
{
  assert(header && truth);

  if(header->detector != int(Detector::kFar)) return 1;

  const int from = TMath::Abs(truth->nonOscNuFlavor);
  const int to = TMath::Abs(truth->nuFlavor);

  if(from == 14 && to == 12 && !useNuE) return 0;

  return pars->TransitionProbability(from, to,
                                     NCType::EEventType(truth->interactionType),
                                     NCType::kBaseLineFar,
                                     truth->nuEnergy);
}


//---------------------------------------------------------------------
void NCEventInfo::SetEventWeight(const bool* useParameters,
                                 const vector<double>& adjustedValues,
                                 const NC::OscProb::OscPars* pars,
                                 SKZPWeightCalculator* skzpcalc,
                                 NC::RunUtil::ERunType runType,
                                 bool useNue,
                                 bool useOscProb)
{
  double neugenWeight = 1.;
  double skzpWeight=GetSKZPWeight(runType);

  //set the defaults for the systematic parameters from the config
  //set the adjustment to be an actual value, the macro should have
  //passed in a number of sigma
  bool use = false;
  for(int i = 0; i < NCType::kNumNeugenParameters; ++i){
    if(useParameters[i]){use = true; break;}
  }

  if(use){
    //Only apply MODBYRS for Carrot or earlier.
    /* StringToType is very heavy on the string comparisons...
       const char* mcString = (header->softVersion).Data();
       ReleaseType::Release_t mcType = ReleaseType::StringToType(mcString);
       bool useMOD = ReleaseType::IsCarrot(mcType);
    */
    const bool useMOD = string(header->softVersion).find("Carrot") != string::npos;
    neugenWeight = FindNeugenWeight(useParameters, adjustedValues, useMOD);
  }

  // Check all the remaining parameters
  for(int i = NCType::kNumNeugenParameters; i < NCType::kNumSystematicParameters; ++i){
    if( useParameters[i] ){use = true; break;}
  }//end loop over parameters

  double recoWeight = 1;
  double crossSectionWeight = 1;
  if(use){
    crossSectionWeight = FindCrossSectionWeight(useParameters, adjustedValues);
    // We need the systematically shifted energies before applying the
    // other systematic weights
    ShiftEnergies(useParameters, adjustedValues);
    recoWeight = FindRecoWeight(useParameters, adjustedValues);
  }

  if(useParameters[NCType::kSKZP]){
    skzpWeight *=
      FindMEGAFitWeightUncertainty(beam->beamType,
                                   runType,
                                   adjustedValues[NCType::kSKZP],
                                   skzpcalc ? skzpcalc : GetSKZPCalc());
  }

  const double survivalProb = useOscProb ?
    FindTransitionProbability(pars, useNue) : 1;

  reco->weight = neugenWeight*recoWeight*crossSectionWeight*skzpWeight*survivalProb;
}

//-----------------------------------------------------------------------
double NCEventInfo::GetSKZPWeight(NC::RunUtil::ERunType runType) const
{
  switch(runType){
  case NC::RunUtil::kRunI:
    return reco_nom->weight;
  case NC::RunUtil::kRunII:
    return reco_nom->weightRunII;
  case NC::RunUtil::kRunIII:
    return reco_nom->weightRunIII;
  default:
    assert(0 && "No weight for requested run");
  }
}

//-----------------------------------------------------------------------
bool NCEventInfo::FinalEventCheck(ANtpAnalysisInfo* ana,
                                  ANtpRecoInfo* rec) const
{
  if(!ana) ana = analysis;
  if(!rec) rec = reco_nom;

  assert(ana && header && rec);

  if(ana->isCC > 0 &&
     ana->isNC < 1 &&
     rec->trackMomentum > 0)
    return false;

  if(!header->isGoodData) return false;

  return true;
}


//-----------------------------------------------------------------------
void NCEventInfo::SetChainBranches(TChain* ch,
                                   TString extractionName,
                                   bool setTruthBranch)
{
  const TString extractionBranch = "analysis"+extractionName;

  ch->SetBranchStatus("*", 0);
  ch->SetBranchStatus("header.*", 1);
  ch->SetBranchStatus("beam.*", 1);
  ch->SetBranchStatus("reco.*", 1);
  ch->SetBranchStatus(extractionBranch+"*", 1);

  ch->SetBranchAddress("header.", &header);
  ch->SetBranchAddress("beam.", &beam);
  // Would like to set this branch to also fill the "reco" pointer,
  // but can't, hence the FillFromChain function
  ch->SetBranchAddress("reco.", &reco_nom);
  ch->SetBranchAddress(extractionBranch, &analysis);

  if(setTruthBranch){
    ch->SetBranchStatus("truth.*", 1);
    ch->SetBranchAddress("truth.", &truth);
  }
}

//-----------------------------------------------------------------------
void NCEventInfo::FillFromChain(TChain* ch, int entry)
{
  ch->GetEntry(entry);
  if(!reco) reco=new ANtpRecoInfo;
  *reco=*reco_nom;
}


//-----------------------------------------------------------------------
void NCEventInfo::SetInfoObjects(ANtpHeaderInfo*    headerInfo,
                                 ANtpBeamInfo*      beamInfo,
                                 ANtpEventInfoNC*   eventInfo,
                                 ANtpTrackInfoNC*   trackInfo,
                                 ANtpShowerInfoNC*  showerInfo,
                                 ANtpAnalysisInfo*  analysisInfo,
                                 ANtpRecoInfo*      recoInfo,
                                 ANtpTruthInfoBeam* truthInfo)
{
  header = headerInfo;
  beam = beamInfo;
  event = eventInfo;
  track = trackInfo;
  shower = showerInfo;
  analysis = analysisInfo;
  reco = recoInfo;
  truth = truthInfo;
}

//-----------------------------------------------------------------------
VldTimeStamp NCEventInfo::GetTimestamp() const
{
  return VldTimeStamp(header->year, header->month, header->day,
                      header->hour, header->minute,
                      TMath::FloorNint(header->second));
}

---