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.24 2009/03/04 12:07:46 bckhouse Exp $");

const ReleaseType::Release_t kReleaseTypeData = ReleaseType::kCedarPhy;
const ReleaseType::Release_t kReleaseTypeMC = ReleaseType::kCedarPhy+ReleaseType::kDaikon;

static const int kQE = 0;
static const int kRES = 1;
static const int kDIS = 2;
static const int kBeamToZbeam[14] = {0, 2, 0, 0, 3, 4, 5, 6, 7, 8, 2, 2, 2, 0};

// 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("NCAnalysisUtils", 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("NCAnalysisUtils", 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];
}

//---------------------------------------------------------------------
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(shower,   ANtpShowerInfoNC);
  COPYFIELD(track,    ANtpTrackInfoNC);
  COPYFIELD(truth,    ANtpTruthInfoBeam);

#undef COPYFIELD
}

//---------------------------------------------------------------------
double NCEventInfo::GetEventVertex(double& x, double& y, double& z)
{
  // 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(reco || (event && track && shower));

  if(reco){
    x = reco->vtxX;
    y = reco->vtxY;
    z = reco->vtxZ;
  }
  else{
    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;
    }
  }

  if(header->detector == int(Detector::kNear)){
    x -= NCType::kNDBeamCenterX;
    y -= NCType::kNDBeamCenterY + NCType::kNDBeamAngle*z;
  }

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


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

  if(reco){
    return reco->nuEnergy;
  }
  else{
    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)
{
  assert(header);

  ReleaseType::Release_t releaseType = kReleaseTypeData;
  SimFlag::SimFlag_t sim = SimFlag::kData;
  if(header->dataType == int(SimFlag::kMC)){
    sim = SimFlag::kMC;
    releaseType = kReleaseTypeMC;
  }

  Detector::Detector_t det = Detector::kNear;
  if(header->detector == int(Detector::kFar))
    det = Detector::kFar;

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

  bool fromCurve = false;
  bool fromRange = false;

  double inputMomentum;

  //check if this is a stopping muon first
  if(reco){
    if(reco->isFullyContained){
      fromRange = true;
      inputMomentum = reco->trackRange;
    }
    else{
      fromCurve = true;
      inputMomentum = reco->trackMomentum;
    }
  }
  // check if this is a stopping muon first or the fit value is stupid,
  // use the range momentum
  else 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, releaseType);
  if(fromRange)
    trackMomentum = FullyCorrectMomentumFromRange(inputMomentum,
                                                  vldc, releaseType);

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


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

  ReleaseType::Release_t releaseType = kReleaseTypeData;
  SimFlag::SimFlag_t sim = SimFlag::kData;
  if(header->dataType == int(SimFlag::kMC)){
    sim = SimFlag::kMC;
    releaseType = kReleaseTypeMC;
  }

  Detector::Detector_t det = Detector::kNear;
  if(header->detector == int(Detector::kFar))
    det = Detector::kFar;

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

  if(reco){
    if(reco->showerEnergy < 0) return ANtpDefaultValue::kDouble;

    return FullyCorrectShowerEnergy(reco->showerEnergy,
                                    showerType, vldc, releaseType);
  }

  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)
    return FullyCorrectShowerEnergy(showerEnergy,
                                    showerType, vldc, releaseType);

  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("NCAnalysisUtils", 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("NCAnalysisUtils", Msg::kWarning)
      << "mcreweight IsNaN "
      <<  truth->nuEnergy
      << " " << truth->interactionType
      << " " << truth->resonanceCode
      << " " << truth->nuFlavor
      << " " << truth->initialState
      << " " << truth->hadronicFinalState
      << " " << nucleus
      << endl;
    weight = 1.;
  }

  if(weight != 1.){
    MAXMSG("NCAnalysisUtils", 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(SKZPWeightCalculator* SKZPWeight, int beamType, int 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;
  if(beam) beamt = beam->beamType;
  SKZPWeightCalculator::RunPeriod_t rp = SKZPWeightCalculator::kRunI;
  if(runType == NC::RunUtil::kRunII) rp = SKZPWeightCalculator::kRunII;

  double pt = TMath::Sqrt(SQR(truth->targetParentPX)
                          +SQR(truth->targetParentPY));
  double emuNew = 0.;
  double eshwNew = 0.;
  double beamweight = 0.;
  double detweight = 0.;
  double weightbeam = SKZPWeight->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("NCAnalysisUtils", Msg::kDebug, 10)
    << "MODBYRS3 = " << weightMODBYRS
    << " beam = " << weightbeam
    << endl;

  return weightMODBYRS*weightbeam;
}


//---------------------------------------------------------------------------
double NCEventInfo::FindMEGAFitWeightUncertainty(int beamType,
                                                 int 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,
                                truth->nuFlavor,
                                truth->nuEnergy,
                                SKZPWeightCalculator::kTotalErrorAfterTune,
                                sigma);
}


//---------------------------------------------------------------------------
double NCEventInfo::FindRecoWeight(const bool* useParameters,
                                   const vector<double>& adjustedValues,
                                   double& energy, double& trackEnergy,
                                   double& showerEnergy) const
{
  assert(analysis);

  //  assert(useParameters.size() == adjustedValues.size());


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

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

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

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

    if(useParameters[NCType::kRelativeHadronicCalibration] &&
       header->detector == Detector::kFar)
      showerEnergy  *= 1. + adjustedValues[NCType::kRelativeHadronicCalibration];

    assert(showerEnergy >= 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)
    showerEnergy  = SimulateShowerOffset(sfINukeHists.offset[offset],
                                         adjustedValues[NCType::kShowerEnergyOffset]);

  if(analysis->isNC>0) energy = showerEnergy;
  if(analysis->isCC>0) energy = trackEnergy + showerEnergy;

  double weight = 1;

  //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.;

  if(useParameters[NCType::kLowCompleteness]
     && analysis->isNC > 0
     && truth->eventCompleteness < 0.5) {
    weight *= 1.+adjustedValues[NCType::kLowCompleteness];
    //cout << "LC w = " << weight << endl;
  }

  //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(showerEnergy <= 0.5)
      weight  *= 1. + 0.078*adjustedValues[NCType::kNCFarCleanNoise];
    else if (showerEnergy > 0.5 && showerEnergy <=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(-showerEnergy)/6.96)*adjustedValues[NCType::kNCFarCleanCR];
  }


  //Tobi's NCNearClean +-15.2% if     Ereco<0.5 GeV
  //                   +- 2.9% if 0.5<EReco<1.0 GeV
  //                   +- 0.4% if 1.0<EReco<1.5 GeV
  if(useParameters[NCType::kNCNearClean]
     && header->detector==Detector::kNear
     && analysis->isNC > 0){
    if(showerEnergy <= 0.5)
      weight  *= 1. + 0.152*adjustedValues[NCType::kNCNearClean];
    else if (showerEnergy > 0.5 && showerEnergy <=1.0)
      weight  *= 1. + 0.029*adjustedValues[NCType::kNCNearClean];
    else if (showerEnergy > 1.0 && showerEnergy <=1.5)
      weight  *= 1. + 0.004*adjustedValues[NCType::kNCNearClean];
    else weight *= 1.;
  }


  //NCCleanRunDiff survey  +-4% if     Ereco<1 GeV in ND
  if(useParameters[NCType::kNCCleanRunDiff]
     && header->detector==Detector::kNear
     && analysis->isNC > 0){
    if(showerEnergy <= 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(showerEnergy <= 1.0)
      weight  *= 1. + 0.04*adjustedValues[NCType::kNCRunDiff];
    else if(showerEnergy > 1.0 && showerEnergy < 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("NCAnalysisUtils", Msg::kDebug, 1) << adjustedValues[i] << " ";
  }
  MAXMSG("NCAnalysisUtils", Msg::kDebug, 1) << endl;

  MAXMSG("NCAnalysisUtils", 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
{
  assert(truth && reco);

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

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

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

    return reco->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->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->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 xfptWeight   = reco->weight;
  if(runType == NC::RunUtil::kRunII) xfptWeight = reco->weightRunII;

  //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 energy = reco->nuEnergy;
  double trackEnergy = reco->muEnergy;
  double showerEnergy = reco->showerEnergy;
  double showerEnergyNC = reco->showerEnergyNC;
  double showerEnergyCC = reco->showerEnergyCC;

  double recoWeight = 1;
  double crossSectionWeight = 1;
  if(use){
    crossSectionWeight = FindCrossSectionWeight(useParameters, adjustedValues);

    // Only need to shift track energy and total energy once. Other times
    // throw the result away in this variable.
    double junk = 0;

    recoWeight = FindRecoWeight(useParameters, adjustedValues, energy,
                                trackEnergy, showerEnergy);
    // Do this again to adjust the NC shower energy too
    FindRecoWeight(useParameters, adjustedValues, junk,
                   junk, showerEnergyNC);
    // ...and the CC shower energy
    FindRecoWeight(useParameters, adjustedValues, junk,
                   junk, showerEnergyCC);
  }

  reco->nuEnergy = energy;
  reco->muEnergy = trackEnergy;
  reco->showerEnergy = showerEnergy;
  reco->showerEnergyNC = showerEnergyNC;
  reco->showerEnergyCC = showerEnergyCC;
  reco->nuEnergyCC = trackEnergy + showerEnergyCC;
  reco->nuEnergyNC = showerEnergyNC ? showerEnergyNC : ANtpDefaultValue::kDouble;


  if(useParameters[NCType::kSKZP])
    xfptWeight *= FindMEGAFitWeightUncertainty(beam->beamType,
                                               runType,
                                               adjustedValues[NCType::kSKZP],
                                               skzpcalc);

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

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


//-----------------------------------------------------------------------
//this is a catch all for final cuts on whether an event should be used.
bool NCEventInfo::FinalEventCheck() const
{
  assert(analysis && header && reco);

  if(analysis->isCC > 0 &&
     analysis->isNC < 1 &&
     reco->trackMomentum > 0)
    return false;

  if(header->detector == Detector::kNear &&
     !NC::RunUtil::IsNDRunGood(header->run, header->subRun))
    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);
  ch->SetBranchAddress("reco.", &reco);
  ch->SetBranchAddress(extractionBranch, &analysis);

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


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

---