NuCuts.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Dec/2005 onwards
//
// Contact: j.j.hartnell@rl.ac.uk

#include <cmath>

#include "TMath.h"
#include "TH1F.h"
#include "TROOT.h"

#include "Conventions/BeamType.h"
#include "DataUtil/infid.h"
#include "DataUtil/infid_finder.h"
#include "Conventions/Munits.h"
#include "Conventions/ReleaseType.h"

#include "NtupleUtils/NuCuts.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuMCEvent.h"
#include "NtupleUtils/NuLibrary.h"
#include "NtupleUtils/NuXMLConfig.h"

using std::endl;
using std::cout;

CVSID("$Id: NuCuts.cxx,v 1.117 2010/02/07 16:00:33 bckhouse Exp $");

//......................................................................
// Functions to identify NuMuBar PQ and NQ selections
Bool_t NuCuts::IsNMB(const NuEvent &nu) {
  return (nu.anaVersion == NuCuts::kNMB0325Alpha ||
          nu.anaVersion == NuCuts::kNMB0325Bravo ||
          nu.anaVersion == NuCuts::kNMB0325BravoTwo ||
          nu.anaVersion == NuCuts::kNMB0325Charlie ||
          nu.anaVersion == NuCuts::kNMB0325Delta ||
          nu.anaVersion == NuCuts::kNMB0325Echo ||
          nu.anaVersion == NuCuts::kNMB0325ChairSound ||
          nu.anaVersion == NuCuts::kRM1 ||
          nu.anaVersion == NuCuts::kRM2 ||
          nu.anaVersion == NuCuts::kJJE2 ||
          nu.anaVersion == NuCuts::kNMBFree ||
          nu.anaVersion == NuCuts::kRHC);
}

//......................................................................

Bool_t NuCuts::IsNMBPQ(const NuEvent &nu) {
  return IsNMB(nu) && (nu.charge == 1);
}

//......................................................................

Bool_t NuCuts::IsNMBNQ(const NuEvent &nu) {
  return IsNMB(nu) && (nu.charge == -1);
}

//......................................................................

NuCuts::NuCuts()
{
  MSG("NuCuts",Msg::kDebug)
    <<"Running NuCuts Constructor..."<<endl;


  MSG("NuCuts",Msg::kDebug)
    <<"Finished NuCuts Constructor"<<endl;
}

//......................................................................

NuCuts::~NuCuts()
{
  MSG("NuCuts",Msg::kDebug)
    <<"Running NuCuts Destructor..."<<endl;


  MSG("NuCuts",Msg::kDebug)
    <<"Finished NuCuts Destructor"<<endl;
}

//......................................................................

Bool_t NuCuts::IsGoodNumberOfTracks(const NuEvent& nu) const
{
  /*
    if (nu.anaVersion==NuCuts::kJJH1 ||
    nu.anaVersion==NuCuts::kJJE1) {
    MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut on number of tracks>=1"<<endl;
    return nu.ntrk>=1;
    }
    else if (nu.anaVersion==NuCuts::kCC0093Std ||
    nu.anaVersion==NuCuts::kCC0250Std ||
    nu.anaVersion==NuCuts::kCC0325Std) {
    MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut on number of tracks>=1"<<endl;
    return nu.ntrk>=1;
    }
    else {
    MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut on number of tracks>=1"<<endl;
    return nu.ntrk>=1;
    }
  */

  // For now, at least, all cuts make the same cut
  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut on number of tracks>=1"<<endl;
  return nu.ntrk >= 1;
}

//......................................................................

Bool_t NuCuts::IsGoodTrackFitPass(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kJJH1) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on trk.fit.pass"<<endl;
    return nu.trkfitpass==1;
  }
  else if (IsNMBPQ(nu))
    {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on trk.fit.pass for +ve tracks"<<endl;
      return nu.trkfitpass==1;
    }
  else if (nu.anaVersion==NuCuts::kJJE1)
    {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on trk.fit.pass for +ve tracks, "
        << "trk.fit.pass w/ trk reclamation for -ve tracks." <<endl;
      if (nu.charge==1){
        // Not currently using track reclamation for RHC - but we could
        //if (nu.detector==Detector::kNear){
        //  return this->IsGoodTrackFitPassReclamation(nu);
        //}
        //else if (nu.detector==Detector::kFar){
        //  return nu.trkfitpass==1;
        //}
        //else {
        //  cout<<"Ahhhh det not known"<<endl;
        //  return nu.trkfitpass==1;
        //}
        return nu.trkfitpass==1;
      }
      else if (nu.charge==-1)
        {
          if (nu.detector==Detector::kNear){
            return this->IsGoodTrackFitPassReclamation(nu);
          }
          else if (nu.detector==Detector::kFar){
            return nu.trkfitpass==1;
          }
          else {
            cout<<"Ahhhh det not known"<<endl;
            return nu.trkfitpass==1;
          }
        }
      else
        {
          cout<<"Bad charge"<<endl;
          return false;
        }
    }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on trk.fit.pass for +ve tracks, "
      << "letting IsNMBNQ deal with -ve tracks." <<endl;
    return nu.trkfitpass==1;
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test || IsNMBNQ(nu))
    {
      if (nu.detector==Detector::kNear) {
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Cutting on trk.fit.pass w/ trk reclamation"<<endl;
        return this->IsGoodTrackFitPassReclamation(nu);
      }
      else if (nu.detector==Detector::kFar) {
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Cutting on trk.fit.pass"
          <<", NOT doing track reclamation for the FD"<<endl;
        return nu.trkfitpass==1;
      }
      else {
        cout<<"Ahhhh det not known"<<endl;
        return nu.trkfitpass==1;
      }
    }
  else {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on trk.fit.pass"<<endl;
    return nu.trkfitpass==1;
    }
}

//......................................................................

Bool_t NuCuts::IsGoodTrackFitPassReclamation(const NuEvent& nu)
{
  // If the track passes the fitter, it passes
  if (nu.trkfitpass==1) return true;

  // Reclaim only +ve's in RHC and only -ve's in FHC
  if (nu.hornIsReverse != (nu.charge == +1)) {
    return false;
  }

  //got these numbers from Masaki's spectrum validation talks
  Bool_t goodTrack= ( abs(nu.planeTrkBegu-nu.planeTrkBegv)<=5  &&
                      abs(nu.planeTrkEndu-nu.planeTrkEndv)<=40 &&
                      nu.planeTrkEndu<270 &&
                      nu.planeTrkEndv<270 );
  return goodTrack;
}

//......................................................................

Bool_t NuCuts::IsGoodPID(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kJJH1) {
    return this->IsGoodAbID(nu);
  }
  else if (nu.anaVersion==NuCuts::kJJE1){
    if (nu.charge==1){
      return this->IsGoodDpID(nu);
    }
    else if (nu.charge==-1){
      return this->IsGoodAbID(nu);
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    return this->IsGoodDpID(nu);
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std) {
    return this->IsGoodDpID(nu);
  }
  else if (nu.anaVersion==NuCuts::kCC0250Std) {
    return this->IsGoodAbID(nu);
  }
  else if (nu.anaVersion==NuCuts::kCC0325Std
           || nu.anaVersion == NuCuts::kCC0720Std
            || IsNMBNQ(nu)
            || nu.anaVersion == NuCuts::kNMB0325ChairSound
            || nu.anaVersion == NuCuts::kMSRock) {
    return this->IsGoodRoID(nu);
  }
  else if (nu.anaVersion==NuCuts::kCC0720Test){
    return (this->IsGoodRoID(nu) || this->IsGoodJmID(nu));
  }
  else if (((nu.anaVersion==NuCuts::kNMB0325Alpha) ||
            (nu.anaVersion==NuCuts::kNMB0325Charlie) ||
            (nu.anaVersion==NuCuts::kNMB0325Delta) ||
            (nu.anaVersion==NuCuts::kRM1) ||
            (nu.anaVersion==NuCuts::kRM2) ||
            (nu.anaVersion==NuCuts::kRHC) ) && IsNMBPQ(nu)) {
    // Both nmb and CC selectors are the same here
    return this->IsGoodRoID(nu);
  }
    else if ((nu.anaVersion == NuCuts::kNMB0325Bravo ||
            nu.anaVersion == NuCuts::kNMB0325BravoTwo) && IsNMBPQ(nu)) {
    return this->IsGoodDpID(nu);
    }
  else if ((nu.anaVersion == NuCuts::kNMB0325Echo) && IsNMBPQ(nu)) {
    return this->IsGoodPoKIN(nu);
  }
  else if ((nu.anaVersion == NuCuts::kNMBFree)&& IsNMBPQ(nu)){
    return true;
  }
  else {
    return this->IsGoodAbID(nu);
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodDpID(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kJJH1) {
    Float_t cutValue=0.4;//my numubar cut
    //Float_t cutValue=0.1;//get away from NC disagreement
    //Float_t cutValue=0.0;//my less harsh numubar cut
    //Float_t cutValue=-0.4;//MK's cut

    if (nu.charge==+1){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on +ve track: IsGoodDpID>"<<cutValue<<endl;
      if (nu.dpID>cutValue) return true;
      else return false;
    }
    else if (nu.charge==-1){
      if (nu.detector==Detector::kNear) cutValue=-0.1;
      else if (nu.detector==Detector::kFar) cutValue=-0.2;
      else {
        cout<<"Ahhh, detector="<<nu.detector<<endl;
        return false;
      }

      //override cutValue
      //cutValue=0.1;//get away from NC disagreement

      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on -ve track: IsGoodDpID>"<<cutValue<<endl;
      if (nu.dpID>cutValue) return true;
      else return false;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE1){
    if (nu.charge==1){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making combined dpID and track jitter cut" <<endl;
      if (nu.dpID < 0.0){return false;}
      if (nu.jitter > 0.5){return false;}
      if (nu.jitter > (0.59*nu.dpID + 0.2)){return false;}
      return true;
    }
    else if(nu.charge==-1){
      Float_t cutValue=-0.1;
      if (nu.detector==Detector::kNear) cutValue=-0.1;
      else if (nu.detector==Detector::kFar) cutValue=-0.2;
      else {
        cout<<"Ahhh, detector="<<nu.detector<<endl;
        return false;
      }
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut: IsGoodDpID>"<<cutValue<<endl;

      if (nu.dpID>cutValue) return true;
      else return false;
    }
    else{
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making combined dpID and track jitter cut" <<endl;
    if (nu.dpID < 0.0){return false;}
    if (nu.jitter > 0.5){return false;}
    if (nu.jitter > (0.59*nu.dpID + 0.2)){return false;}
    return true;
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std) {
    Float_t cutValue=-0.1;
    if (nu.detector==Detector::kNear) cutValue=-0.1;
    else if (nu.detector==Detector::kFar) cutValue=-0.2;
    else {
      cout<<"Ahhh, detector="<<nu.detector<<endl;
      return false;
    }
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodDpID>"<<cutValue<<endl;

    if (nu.dpID>cutValue) return true;
    else return false;
  }
  else if (nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std || 
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test || IsNMBNQ(nu)) {
    Float_t cutValue=-0.1;
    if (nu.detector==Detector::kNear) cutValue=-0.1;
    else if (nu.detector==Detector::kFar) cutValue=-0.2;
    else {
      cout<<"Ahhh, detector="<<nu.detector<<endl;
      return false;
    }
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodDpID>"<<cutValue<<endl;

    if (nu.dpID>cutValue) return true;
    else return false;
  }
   else if ((nu.anaVersion == NuCuts::kNMB0325Bravo ||
            nu.anaVersion == NuCuts::kNMB0325BravoTwo) && IsNMBPQ(nu)) {
    Float_t cutValue=0.25;

    if (nu.anaVersion == NuCuts::kNMB0325BravoTwo) {
      cutValue = 0.32;
      }
    MAXMSG("NuCuts",Msg::kInfo,1)
      << "Making +ve  cut: IsGoodDpID>"<<cutValue<<endl;
    if (nu.dpID>cutValue)
      return true;
    else
      return false;
      }
  else {
    //Float_t cutValue=0.4;//my numubar cut
    Float_t cutValue=0.0;//my less harsh numubar cut
    //Float_t cutValue=-0.4;//MK's cut

    if (nu.charge==+1){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on +ve track: IsGoodDpID>"<<cutValue<<endl;
      if (nu.dpID>cutValue) return true;
      else return false;
    }
    else if (nu.charge==-1){
      if (nu.detector==Detector::kNear) cutValue=-0.1;
      else if (nu.detector==Detector::kFar) cutValue=-0.2;
      else {
        cout<<"Ahhh, detector="<<nu.detector<<endl;
        return false;
      }
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on -ve track: IsGoodDpID>"<<cutValue<<endl;

      if (nu.dpID>cutValue) return true;
      else return false;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodRoID(const NuEvent& nu) const
{
  Float_t cutValue=0.3;

  // Check for the NuMuBar cuts first, then drop down to the default (original)
  // routine if none of them are to be used.
  if (IsNMBPQ(nu)) {
    if ((nu.anaVersion == NuCuts::kNMB0325Alpha) ||
        (nu.anaVersion == NuCuts::kNMB0325Charlie) ||
        (nu.anaVersion == NuCuts::kRM1)) {
      cutValue = 0.625;
    }
    else if (nu.anaVersion == NuCuts::kRM2) {
      cutValue = 0.80;
    }
    else if (nu.anaVersion == NuCuts::kNMB0325Delta) {
      cutValue = 0.55;
    }
    else if (nu.anaVersion == NuCuts::kRHC) {
        cutValue = 0.30; // Prelim. RHC cut the same as the CC cut
    }
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making +ve IsGoodRoID>"<<cutValue<<endl;
  }

  // Drop to original behaviour in two cases:
  // - Is NOT a numubar selector
  // - IS a nmb selector, but this event has negative charge (i.e. a neutrino event)
  if (!IsNMB(nu) || IsNMBNQ(nu)) {
    // Original Behaviour
    if (nu.detector==Detector::kNear) cutValue=0.3;
    else if (nu.detector==Detector::kFar) cutValue=0.3;
    else {
      cout<<"Ahhh, detector="<<nu.detector<<endl;
      return false;
    }
  }

  // Rock muons are different!
  if (nu.anaVersion == NuCuts::kMSRock) {
    cutValue = 0.20;
  }

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut: IsGoodRoID>"<<cutValue<<endl;

  return (nu.roID>cutValue);
  return true;
}
Bool_t NuCuts::IsGoodJmID(const NuEvent& nu) const
{
  Float_t cutValue=0.5;
  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut: IsGoodJmID>"<<cutValue<<endl;
  return (nu.jmID>cutValue);
  return true;
} 
Bool_t NuCuts::IsGoodNtID(const NuEvent& /*nu*/) const
{
  Float_t cutValue=0.8;
  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making cut: ISGoodNtID>"<<cutValue<<endl;
  //  return (nu.ntID>cutValue);
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodPoKIN(const NuEvent& nu) const
{
  Float_t cutValue=-0.1;
  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Making Cut: IsGoodPoKIN>" << cutValue << endl;
  // Check the value!
  if (nu.anaVersion == NuCuts::kNMB0325Echo) {
    if (nu.poIDKin <= cutValue) return false;
  }

  // return true for all other cases (not cutting, is correct)
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodAbID(const NuEvent& nu) const
{
  Float_t cutValue=0.85;

  if (nu.anaVersion==NuCuts::kJJH1) {
    cutValue=0.99;

    if (nu.charge==+1){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on +ve track: IsGoodAbID>"<<cutValue<<endl;
      if (nu.abID>cutValue) return true;
      else return false;
    }
    else if (nu.charge==-1){
      if (nu.detector==Detector::kNear) cutValue=0.85;
      else if (nu.detector==Detector::kFar) cutValue=0.85;
      else {
        cout<<"Ahhh, detector="<<nu.detector<<endl;
        return false;
      }

      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on -ve track: IsGoodAbID>"<<cutValue<<endl;
      if (nu.abID>cutValue) return true;
      else return false;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE1){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodAbID>"<<cutValue<<endl;
    if (nu.abID>cutValue) return true;
    else return false;
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodAbID>"<<cutValue<<endl;
    if (nu.abID>cutValue) return true;
    else return false;
  }
  else if (nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test ||  IsNMBNQ(nu)) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodAbID>"<<cutValue<<endl;
    if (nu.abID>cutValue) return true;
    else return false;
  }
  else {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut: IsGoodAbID>"<<cutValue<<endl;
    if (nu.abID>cutValue) return true;
    else return false;
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodSigmaQP_QP(const NuEvent& nu) const
{
  Float_t cutValue=0.3;

  if (nu.anaVersion==NuCuts::kJJH1) {
    //first cut on anti-neutrinos
    if (nu.charge>0){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on +ve track: SigmaQP/QP<"<<cutValue<<endl;
      if (fabs(nu.sigqp_qp)<cutValue) {
        return true;
      }
      else return false;
    }
    //now look at neutrinos
    else if (nu.charge<0) {
      Bool_t doCut=false;//switch
      if (doCut) {
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"If curv used, making cut on -ve track: SigmaQP/QP>"
          <<-1*cutValue<<"; if range used, no cut"<<endl;
        if (nu.usedRange) {
          return true;//don't cut on fit error if using range
          //neutrinos swamp anti-neutrinos so little background
        }
        else if (nu.usedCurv) {
          if (fabs(nu.sigqp_qp)<cutValue) return true;
          else return false;
        }
        else {
          MAXMSG("NuCuts",Msg::kWarning,10)
            <<"Ahh, usedRange="<<nu.usedRange
            <<"and usedCurv="<<nu.usedCurv<<endl;
          return false;
        }
      }
      else {
        //don't cut
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Not cutting on SigmaQP/QP for -ve tracks"<<endl;
        return true;
      }
    }
    else {
      cout<<"Ahh, bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE1){
    if (1==nu.charge){
      cutValue = 0.065;
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on sigmaQP>" << cutValue <<endl;
      if (nu.sigqp>cutValue){return false;}
      else{return true;}
    }
    else if (-1==nu.charge){
      //no cut
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Not cutting on SigmaQP/QP"<<endl;
      return true;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    cutValue = 0.065;
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut on sigmaQP>" << cutValue <<endl;
    if (nu.sigqp>cutValue){return false;}
    else{return true;}
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test || IsNMBNQ(nu)) {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on SigmaQP/QP"<<endl;
    return true;
  }
    else if (((nu.anaVersion == NuCuts::kNMB0325Bravo) ||
            (nu.anaVersion == NuCuts::kNMB0325Echo)) && IsNMBPQ(nu)){
    cutValue = 3.5;
    MAXMSG("NuCuts",Msg::kInfo,1)
      << "Making Cut on +ve, QP_sigmaQP>" << cutValue << endl;
    if (1./nu.sigqp_qp > cutValue)
      return true;
    else
      return false;

      }
  else if ((nu.anaVersion == NuCuts::kRM2) && IsNMBPQ(nu)){
    cutValue = 2.3;
    MAXMSG("NuCuts",Msg::kInfo,1)
      << "Making Cut on +ve, QP_sigmaQP>" << cutValue << endl;
    if (1./nu.sigqp_qp > cutValue)
      return true;
    else
      return false;
  }
  else if ((nu.anaVersion == NuCuts::kRM1) && IsNMBPQ(nu)) {
    cutValue = 2.0;
    MAXMSG("NuCuts",Msg::kInfo,1)
      << "Making Cut on +ve, QP_sigmaQP>" << cutValue << endl;
    if (1./nu.sigqp_qp > cutValue)
      return true;
    else
      return false;
  }
  else {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on SigmaQP/QP"<<endl;
    return true;
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodRelAngle(const NuEvent &nu) const
{
  // Calculate relative angle - pi
  Float_t cutValue = 0;

  // Only cut on NMB Selectors
  if (!IsNMB(nu))
    return true;

  // Don't cut on NMB negative charge events
  if (nu.charge == -1)
    return true;
  else if (nu.charge != 1)
    {
      cout << "Bad Charge: " << nu.charge << endl;
      return false;
    }

  // Do Rashid's relA selector for RM1
  if (nu.anaVersion == NuCuts::kRM1)
    {
      if (nu.relativeAngle < 1.0 || nu.relativeAngle > 5.0)
        return true;
      else
        return false;
    }

  switch (nu.anaVersion)
    {
    case NuCuts::kNMB0325Alpha:
      cutValue = 2.12;
      break;
       case NuCuts::kNMB0325Bravo:
      cutValue = 2.12;
      break;
    case NuCuts::kNMB0325BravoTwo:
      cutValue = 2.12;
      break;
    case NuCuts::kNMB0325Echo:
      cutValue = 2.02;
      break;
    case NuCuts::kRM2:
      cutValue = 2.0;
      break;
    default:
      // Do not cut if not one of the above
      MAXMSG("NuCuts",Msg::kInfo,1)
        << "Not Cutting on |relA-Pi|" << endl;
      return true;
    }

  MAXMSG("NuCuts",Msg::kInfo,1)
    << "Cutting on |relA-Pi|>" << cutValue << endl;
  if (abs(nu.relativeAngle - TMath::Pi()) <= cutValue)
    return false;

  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodTrackLength(const NuEvent &nu) const
{
  Int_t cutValue = 0;

  // Only cut on NMB Selectors
  if (!IsNMB(nu))
    return true;

  // Don't cut on NMB negative charge events
  if (nu.charge == -1)
    return true;
  else if (nu.charge != 1)
    {
      cout << "Bad Charge: " << nu.charge << endl;
      return false;
    }

  switch (nu.anaVersion)
    {
    case NuCuts::kNMB0325Alpha:
    case NuCuts::kNMB0325Charlie:
      cutValue = 35;
      break;
    case NuCuts::kNMB0325Delta:
      cutValue = 37;
      break;
    case NuCuts::kRM1:
      cutValue = 36;
      break;
    case NuCuts::kJJE2:
      cutValue = 34;
      break;
    default:
      // Do not cut if not one of the above
      MAXMSG("NuCuts",Msg::kInfo,1) << "Not Cutting on trkL" << endl;
      return true;
    }

  // Since this function only cuts for NuMuBar selectors at the moment,
  // we can make a blanket acceptance for NQ events
  if (nu.charge == -1)
    return true;
  else if (nu.charge != 1) {
    cout << "Bad Charge: " << nu.charge << endl;
    return false;
  }

  MAXMSG("NuCuts",Msg::kInfo,1)
    << "Cutting on trkL>" << cutValue << endl;
  if (nu.trkLength <= cutValue)
    return false;

  return true;
}

//......................................................................


Bool_t NuCuts::IsGoodTrackLengthInUV(const NuEvent &nu) const
{
  Int_t cutValue = 19;
  Int_t DiffU=0;
  Int_t DiffV=0;

  // Only cut on NMB Selectors
  if (!IsNMB(nu))
    return true;

  // Don't cut on NMB negative charge events
  if (nu.charge == -1)
    return true;
  else if (nu.charge != 1)
    {
      cout << "Bad Charge for TrkLenghInUV cut: " << nu.charge << endl;
      return false;
    }

  if (nu.anaVersion == NuCuts::kRM2) {
    DiffV=nu.planeTrkEndv -nu.planeTrkBegv;
    DiffU=nu.planeTrkEndu -nu.planeTrkBegu;
    MAXMSG("NuCuts",Msg::kInfo,1)
      << "Making Cut on UV Track length >" << cutValue << endl;
    if (DiffU > cutValue && DiffV > cutValue)
      return true;
    else
      return false;
  }
  else {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on UV Track length"<<endl;
    return true;
  }
}

//......................................................................

Bool_t NuCuts::IsGoodQP(const NuEvent& nu) const
{
  MAXMSG("NuCuts",Msg::kWarning,1) << "CUTTING ON QP. YOU PROBABLY DON'T WANT TO DO THIS" << endl;

  switch (nu.anaVersion)
    {
    case NuCuts::kNMB0325Alpha:
    case NuCuts::kNMB0325Bravo:
    case NuCuts::kNMB0325BravoTwo:
    case NuCuts::kNMB0325Charlie:
    case NuCuts::kNMB0325Delta:
    case NuCuts::kNMB0325Echo:
    case NuCuts::kRM2:
    case NuCuts::kNMBFree:
    case NuCuts::kRHC:
      { // brackets defeat "jump to case crosses intialization" errors
        MAXMSG("NuCuts",Msg::kInfo,1) << "Cutting on Q/P>0" << endl;
        if (nu.qp <= 0)
          return false;
        break;
      }
    default:
      {
        MAXMSG("NuCuts",Msg::kInfo,1) << "Not Cutting on Q/P" << endl;
      }
    }

  // Either not cutting or passed cut
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodFitChi2PerNdof(const NuEvent& nu) const
{
  Float_t cutValue=10;

  if (nu.anaVersion==NuCuts::kJJH1) {
    Bool_t doCut=false;
    if (doCut) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on chi2PerNdof<"<<cutValue<<endl;
      //make cut
      if (nu.chi2PerNdof>cutValue) {
        return false;
      }
      else return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Not cutting on Chi2PerNdof"<<endl;
      return true;
    }
  }

  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kJJE1 ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test || IsNMB(nu)) {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Chi2PerNdof"<<endl;
    return true;
  }
  else {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Chi2PerNdof"<<endl;
    return true;
  }
}

//......................................................................

Bool_t NuCuts::IsGoodFitProb(const NuEvent& nu) const
{
  Float_t cutValue=0.1;

  if (nu.anaVersion==NuCuts::kJJH1) {
    if (nu.charge==+1){
      cutValue=0.1;
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on +ve track: Fit Prob>"<<cutValue<<endl;
      if (nu.prob>cutValue) return true;
      else return false;
    }
    else if (nu.charge==-1){
      Bool_t doCut=false;

      if (doCut) {
        //if (nu.detector==Detector::kNear) cutValue=-0.1;
        //else if (nu.detector==Detector::kFar) cutValue=-0.2;
        //else {
        //cout<<"Ahhh, detector="<<nu.detector<<endl;
        //return false;
        //}
        cutValue=0.0;

        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Making cut on -ve track: Fit Prob>="<<cutValue<<endl;
        if (nu.prob>=cutValue) return true;
        else return false;
      }
      else {
        //no cut
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Not cutting on trk fit probability for -ve tracks"<<endl;
        return true;
      }
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE1){
    if (1==nu.charge){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on tfProb>" << cutValue <<endl;
      cutValue = 0.1;
      if (nu.prob<0.1){return false;}
      else{return true;}
    }
    else if (-1==nu.charge){
      //no cut
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Not cutting on Fit Prob"<<endl;
      return true;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut on tfProb>" << cutValue <<endl;
    cutValue = 0.1;
    if (nu.prob<0.1){return false;}
    else{return true;}
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test ||

           IsNMBNQ(nu) ) {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Fit Prob"<<endl;
    return true;
  }
  else if ((nu.anaVersion == NuCuts::kNMB0325Delta) && IsNMBPQ(nu))
    {
      MAXMSG("NuCuts",Msg::kInfo,1) << "Cutting on prob>0.003" << endl;
      if (nu.prob <= 0.003)
        return false;
      else
        return true;
    } else {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Fit Prob"<<endl;
    return true;
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodMajorityCurvature(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kJJE1){
    if (1==nu.charge){
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Making cut on Majority Curvature < 0" << endl;
      if (nu.smoothMajC<0.0){return false;}
      else{return true;}
    }
    else if (-1==nu.charge){
      //no cut
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Not cutting on Majority Curvature for charge=-1"<<endl;
      return true;
    }
    else {
      cout<<"Bad charge"<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kJJE2 && nu.charge==1){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut on Majority Curvature < 0" << endl;
    if (nu.smoothMajC<0.0){return false;}
    else{return true;}
  }
  else if (nu.anaVersion==NuCuts::kJJH1 ||
           nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test ||
           IsNMBNQ(nu) ) {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Majority Curvature"<<endl;
    return true;
  }
  else if ((nu.anaVersion == NuCuts::kNMB0325Charlie) && IsNMBPQ(nu)) {
    MAXMSG("NuCuts",Msg::kInfo,1)<<"Making cut on Majority Curvature > -0.22" << endl;

    // Make the numubar cut
    if (nu.smoothMajC > -0.22)
      return true;
    else
      return false;
  }
  else {
    //no cut
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on Majority Curvature"<<endl;
    return true;
    }
  return true;
}

//......................................................................

Bool_t NuCuts::IsLI(const NuEvent& nu) const
{
  //for the ND just cut on litime
  if (nu.detector==Detector::kNear) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on ND LI events using litime!=-1"<<endl;
    return (nu.litime!=-1);
  }

  //for the FD cut depending on the analysis version
  if (nu.anaVersion==NuCuts::kJJH1 ||
      nu.anaVersion==NuCuts::kJJE1 ||
      nu.anaVersion==NuCuts::kCC0325Std ||
      nu.anaVersion==NuCuts::kCC0720Std ||
      nu.anaVersion==NuCuts::kCC0720Test ||
      IsNMB(nu)) {
    //cut on LI
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on LI events using LISieve and litime!=-1"<<endl;
    return (nu.isLI || nu.litime!=-1);
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on LI events using FDRCBoundary and litime!=-1"<<endl;
    return (nu.litag>0 || nu.litime!=-1);
  }
  else {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on LI events using LISieve and litime!=-1"<<endl;
    return (nu.isLI || nu.litime!=-1);
  }
}

//......................................................................

Bool_t NuCuts::IsCoilOk(const NuEvent& nu) const
{
  //check if MC or data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on IsCoilOk for simFlag="<<nu.simFlag<<endl;
    return true;
  }

  //check if the flag is set to cut
  if (!nu.cutOnDataQuality) {//use data quality flag
    MAXMSG("NuCuts",Msg::kWarning,1)
      <<"Not cutting on IsCoilOk (flag set to not cut)"<<endl;
    return true;
  }

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on coilIsOk"<<endl;

  //check if coil is ok (can be forward or reversed)
  if (nu.coilIsOk) {
    return true;
  }
  else {
    MAXMSG("NuCuts",Msg::kInfo,10)
      <<"Bad magnetic coil, coilIsOk="<<nu.coilIsOk<<endl;
    return false;
  }
}

//......................................................................

Bool_t NuCuts::IsGoodCoilCurrent(const NuEvent& nu) const
{
  //check if MC or data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on CoilCurrent for simFlag="<<nu.simFlag<<endl;
    return true;
  }

  //check if Near or not
  if (nu.detector!=Detector::kNear) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on ND CoilCurrent for detector="<<nu.detector
      <<endl;
    return true;
  }

  //check if the flag is set to cut
  if (!nu.cutOnDataQuality) {//use data quality flag
    MAXMSG("NuCuts",Msg::kWarning,1)
      <<"Not cutting on CoilCurrent (flag set to not cut)"<<endl;
    return true;
  }

  Float_t valueLow=-1000;
  Float_t valueHigh=1000;

  if (nu.anaVersion==NuCuts::kJJH1 ||
      nu.anaVersion==NuCuts::kJJE1 ||
      IsNMB(nu)) {
    //allow reversed field currents
    if (nu.coilCurrent<valueLow || nu.coilCurrent>valueHigh) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on CoilCurrent < "<<valueLow
        <<" || CoilCurrent > "<<valueHigh<<endl;
      return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,10)
        <<"Bad coil current="<<nu.coilCurrent<<endl;
      return false;
    }
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test) {
    if (nu.coilCurrent<valueLow) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on CoilCurrent < "<<valueLow<<endl;
      return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,10)
        <<"Bad coil current="<<nu.coilCurrent<<endl;
      return false;
    }
  }
  else {
    if (nu.coilCurrent<valueLow) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on CoilCurrent < "<<valueLow<<endl;
      return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,10)
        <<"Bad coil current="<<nu.coilCurrent<<endl;
      return false;
    }
  }
}

//......................................................................

Bool_t NuCuts::IsGoodBeamDetPOTCountingStage(const NuEvent& nu) const
{

  //return if MC
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on GoodBeam for simFlag="<<nu.simFlag<<endl;
    return true;
  }

  //don't do this here, since this is the sntp beam data
  //check if the flag is set to cut
  //if (!nu.cutOnBeamInfo) {
  //MAXMSG("NuCuts",Msg::kWarning,1)
  //    <<"Not cutting on BeamInfo (flag set to not cut)"<<endl;
  //return true;
  //}

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on IsGoodBeamDetPOTCountingStage"<<endl;

  if (nu.goodBeamSntp && nu.coilIsOk) {
    return true;
  }
  return false;
}

//......................................................................

Bool_t NuCuts::IsGoodBeam(const NuEvent& nu) const
{

  //return if MC
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on GoodBeam for simFlag="<<nu.simFlag<<endl;
    return true;
  }

  //check if the flag is set to cut
  if (!nu.cutOnBeamInfo) {
    MAXMSG("NuCuts",Msg::kWarning,1)
      <<"Not cutting on BeamInfo (flag set to not cut)"<<endl;
    return true;
  }

  //Trish's cuts:
  //goodbeam==1&&hornI<-10&&(beamcnf==3||beamcnf==10||beamcnf==11)

  //3=kL010z185i - LE-std
  //7=kL250z200i - phe (removed for CC0250 results)
  //10=kL010z170i - LE lower horn cur
  //11=kL010z200i - LE higher horn cur
  //15=kL150z200i - pmhe (removed for CC0250 results)

  if (nu.anaVersion==NuCuts::kCC0093Std ||
      nu.anaVersion==NuCuts::kCC0250Std) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on GoodBeam with CC0250Std cuts"<<endl;
    if (nu.goodBeam &&
        nu.hornCur<-10 && nu.hornCur!=-999999 &&
        (nu.beamTypeDB==BeamType::kL010z185i ||
         nu.beamTypeDB==BeamType::kL010z170i ||
         nu.beamTypeDB==BeamType::kL010z200i)){
      return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,100)
        <<"Bad beam: goodBeam="<<nu.goodBeam
        <<", goodBeamSntp="<<nu.goodBeamSntp
        <<", hornCur="<<nu.hornCur
        <<", nu.beamTypeDB="<<nu.beamTypeDB
        <<", hornIsReverse="<< (nu.hornIsReverse?"Yes":"No")
        << endl;
      if (nu.detector==Detector::kNear) {
        MAXMSG("NuCuts",Msg::kError,1000)
            <<"IsGoodBeam: cutting on Bad beam here messes"
          <<" up POT counting.  Event: " << nu.run << "/" << nu.subRun << ", " <<  nu.snarl <<endl;
      }
      return false;
    }
  }
  else if (nu.anaVersion == NuCuts::kCC0325Std ||
           nu.anaVersion == NuCuts::kCC0720Std ||
           nu.anaVersion == NuCuts::kCC0720Test || IsNMB(nu)){
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on GoodBeam with CC0325Std cuts"<<endl;


    bool isGood = true;

    // Check current
    if (!nu.goodBeam)
      isGood = false;
    if (nu.hornCur==-999999)
      isGood = false;

    if (nu.hornCur < -10) {
      if (nu.beamTypeDB != BeamType::kL010z185i &&
          nu.beamTypeDB != BeamType::kL010z170i &&
          nu.beamTypeDB != BeamType::kL010z200i &&
          nu.beamTypeDB != BeamType::kL250z200i) {
        isGood = false;
      }
      if (nu.hornIsReverse) {
        MAXMSG("NuCuts",Msg::kError,1000)
        << "nu.hornIsReverse is true but it shouldn't be for beamTypeDB="
        << nu.beamTypeDB << endl;
      }
    }
    else if (nu.hornCur >  10) {
      if (nu.beamTypeDB != BeamType::kL010z185i_rev)
        isGood = false;
      if (!nu.hornIsReverse) {
        MAXMSG("NuCuts",Msg::kError,1000)
        << "nu.hornIsReverse is false but it shouldn't be for beamTypeDB="
        << nu.beamTypeDB << endl;
      }
    }
    else {
      if (nu.beamTypeDB != BeamType::kL010z000i)
        isGood = false;
    }


    if (!isGood) {
      MAXMSG("NuCuts",Msg::kInfo,100)
      <<"Bad beam: goodBeam="<<nu.goodBeam
      <<", goodBeamSntp="<<nu.goodBeamSntp
      <<", hornCur="<<nu.hornCur
      <<", nu.beamTypeDB="<<nu.beamTypeDB
      <<", hornIsReverse="<< (nu.hornIsReverse?"Yes":"No")
      << endl;
      if (nu.detector==Detector::kNear) {
        MAXMSG("NuCuts",Msg::kError,1000)
        <<"IsGoodBeam: cutting on Bad beam here messes"
        <<" up POT counting.  Event: " << nu.run << "/" << nu.subRun << ", " <<  nu.snarl <<endl;
      }
      return false;
    }

    //Make sure we're not getting mixed beam
    static Bool_t firstGoodEvent = true;
    static Int_t firstType = 0;

    if (firstGoodEvent){
      firstGoodEvent = false;
      firstType = nu.beamTypeDB;
      // Treat le10 170, 185, 200 as all the same
      if (firstType == BeamType::kL010z170i ||
          firstType == BeamType::kL010z200i) {
        firstType = BeamType::kL010z185i;
      }
    }
    else{
      Int_t thisType = nu.beamTypeDB;
      // Treat le10 170, 185, 200 as all the same
      if (nu.simFlag==SimFlag::kData && nu.detector==Detector::kFar){
        if (thisType == BeamType::kL010z170i ||
            thisType == BeamType::kL010z200i) {
          thisType = BeamType::kL010z185i;
        }
      }

      if (thisType != firstType) {
        MAXMSG("NuCuts",Msg::kWarning,100)
        << "Seeing mixed beam: firstType=" << firstType
        << ", thisType=" << thisType
        << ", beamTypeDB=" << nu.beamTypeDB << endl;
        if (nu.detector==Detector::kNear) {
          MAXMSG("NuCuts",Msg::kError,1000)
          <<"IsGoodBeam: cutting on Bad beam here messes"
          <<" up POT counting.  Event: " << nu.run << "/" << nu.subRun << ", " << nu.snarl <<endl;
        }
        return false;
      }
    }
    //Beam consistent with first event in DST, so we can let it pass
    return true;
  }
  else {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on GoodBeam"<<endl;
    return true;
  }
}

//......................................................................

Bool_t NuCuts::IsGoodDataQuality(const NuEvent& nu) const
{
  //return if MC
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on DataQuality for simFlag="<<nu.simFlag<<endl;
    return true;
  }

  //check if the flag is set to cut
  if (!nu.cutOnDataQuality) {
    MAXMSG("NuCuts",Msg::kWarning,1)
      <<"Not cutting on DataQuality (flag set to not cut)"
      <<", isGoodDataQuality="<<nu.isGoodDataQuality<<endl;
    return true;
  }

  //make the cut
  if (nu.detector==Detector::kFar) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on IsGoodDataQuality"<<endl;
    return nu.isGoodDataQuality;
  }
  else if (nu.detector==Detector::kNear) {
    if (nu.anaVersion==NuCuts::kCC0093Std ||
        nu.anaVersion==NuCuts::kJJH1 ||
        nu.anaVersion==NuCuts::kCC0250Std ||
        nu.anaVersion==NuCuts::kJJE1 ||
        nu.anaVersion==NuCuts::kCC0325Std ||
        nu.anaVersion == NuCuts::kCC0720Std ||
        nu.anaVersion==NuCuts::kCC0720Test ||
        nu.anaVersion==NuCuts::kJJE2 ||
        nu.anaVersion == NuCuts::kNMB0325Alpha ||
        nu.anaVersion == NuCuts::kNMB0325Bravo ||
        nu.anaVersion == NuCuts::kNMB0325BravoTwo ||
        nu.anaVersion == NuCuts::kNMB0325Charlie ||
        nu.anaVersion == NuCuts::kNMB0325Delta ||
        nu.anaVersion == NuCuts::kNMB0325Echo  ||
        nu.anaVersion == NuCuts::kRM1  ||
        nu.anaVersion == NuCuts::kRM2  ||
        nu.anaVersion == NuCuts::kNMBFree) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"NOT cutting on IsGoodDataQuality for anaVersion="<<nu.anaVersion<<endl;
      //don't cut on coil current or data quality here due to
      //POT counting (now changed so that the "good" pot is
      //after data quality and coil current)
      //the data quality information was not available before the 7e20 analyses
      return true;
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on IsGoodDataQuality"<<endl;
      return nu.isGoodDataQuality;
    }
  }
  else {
    cout<<"Ahhh, det="<<nu.detector<<endl;
    return false;
  }
}

//......................................................................

Bool_t NuCuts::IsGoodTimeToNearestSpill(const NuEvent& nu) const
{

  //return if MC
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on TimeToNearestSpill for simFlag="<<nu.simFlag
      <<endl;
    return true;
  }

  //return if ND
  if (nu.detector!=Detector::kFar) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on TimeToNearestSpill for detector="<<nu.detector
      <<endl;
    return true;
  }

  //check if the flag is set to cut
  if (!nu.cutOnSpillTiming) {
    MAXMSG("NuCuts",Msg::kWarning,1)
      <<"Not cutting on SpillTiming (flag set to not cut)"<<endl;
    return true;
  }

  //define the default time cuts
  Double_t lowerTimeBound=-20*(Munits::microsecond);
  Double_t upperTimeBound=+30*(Munits::microsecond);

  if (this->IsNMB(nu) || NuCuts::kJJE1==nu.anaVersion){
    //New and improved spill timing cut.
    lowerTimeBound=-2*(Munits::microsecond);
    upperTimeBound=+12*(Munits::microsecond);
  }

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on TimeToNearestSpill: "
    <<lowerTimeBound/(Munits::microsecond)<<"<t<"
    <<upperTimeBound/(Munits::microsecond)<<" us"<<endl;

  if (nu.timeEvtMin-nu.timeToNearestSpill>lowerTimeBound &&
      nu.timeEvtMin-nu.timeToNearestSpill<upperTimeBound) {
    return true;
  }
  else return false;
}

//......................................................................

Bool_t NuCuts::IsGoodDirCos(const NuEvent& nu) const
{
  Float_t cutValue=0.6;
  //trkvtxdcosz //dirCosNu
  Float_t value=nu.dirCosNu;

  if (nu.detector==Detector::kNear) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on track direction cosine for ND"<<endl;
    return true;
  }
  else if (nu.detector==Detector::kFar) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making cut on track direction cosine>"<<cutValue<<endl;
    if (value>cutValue) return true;
    else return false;
  }
  else {
    cout<<"Ahhh, detector="<<nu.detector<<endl;
    return false;
  }
}

//......................................................................

Bool_t NuCuts::IsInFidVol(Float_t x,Float_t y,Float_t z,
                          Float_t u,Float_t v,
                          Int_t planeTrkVtx,Float_t rTrkVtx,
                          Int_t detector,Int_t anaVersion,
                          Int_t releaseType,Int_t simFlag) const
{

  static Bool_t infid_initialised=false;

  //NEAR detector
  if (detector==Detector::kNear){
    if (anaVersion==NuCuts::kCC0093Std) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on ND xyz point being in kCC0093Std fid vol"<<endl;
      return this->IsInFidVolNDCC0093Std(x,y,z);
    }
    else if (anaVersion==NuCuts::kCC0250Std) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on ND xyz point being in kCC0250Std fid vol"<<endl;
      return this->IsInFidVolNDCC0250Std(x,y,z);
    }
    else if ( anaVersion==NuCuts::kCC0325Std ||
              anaVersion == NuCuts::kCC0720Std ||
              anaVersion==NuCuts::kCC0720Test ||
              anaVersion==NuCuts::kJJH1 ||
              anaVersion==NuCuts::kJJE1 ||
              anaVersion==NuCuts::kJJE2 ||
              anaVersion == NuCuts::kNMB0325Alpha ||
              anaVersion == NuCuts::kNMB0325Bravo ||
              anaVersion == NuCuts::kNMB0325BravoTwo ||
              anaVersion == NuCuts::kNMB0325Charlie ||
              anaVersion == NuCuts::kNMB0325Delta ||
              anaVersion == NuCuts::kNMB0325Echo  ||
              anaVersion == NuCuts::kRM1  ||
              anaVersion == NuCuts::kRM2  ||
              anaVersion == NuCuts::kNMBFree ||
              anaVersion == NuCuts::kRHC ||
              anaVersion == NuCuts::kMSRock ) {
      static Int_t msgCounter = 0;
      if (msgCounter < 20){
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Cutting on ND xyz point being in kCC0325Std fid vol"<<endl;
        ++msgCounter;
      }
      if (!infid_initialised) {
        infid_initialised=true;
        choose_infid_set("cc2008");
      }
      return infid(static_cast<Detector::Detector_t>(detector),
        static_cast<SimFlag::SimFlag_t>(simFlag),
        x,y,z);
    }
    else if (anaVersion==NuCuts::kFullDST) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"For kFullDST cutting on ND xyz point being in"
        <<" \"loose\" fid vol"<<endl;
      return this->IsInFidVolLoose(x,y,z,detector);
    }
    else {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on ND xyz point being in Pitt fid vol"<<endl;
      return this->IsInFidVolPitt(x,y,z,u,v);
    }
  }
  //FAR detector
  else if (detector==Detector::kFar){
    if (anaVersion==NuCuts::kCC0093Std) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on FD xyz point being in kCC0093Std fid vol"<<endl;
      return this->IsInFidVolFDCC0093Std(x,y,z);
    }
    else if (anaVersion==NuCuts::kCC0250Std) {
      if (ReleaseType::IsBirch(releaseType)) {
        MSG("NuCuts",Msg::kError)
          <<"You can't run kCC0250Std cuts on Birch data since not all"
          <<" the variables are available!"<<endl;
        MSG("NuCuts",Msg::kFatal)<<"This doesn't print"<<endl;
      }
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on FD xyz point being in kCC0250Std fid vol"<<endl;
      return this->IsInFidVolFDCC0250Std(planeTrkVtx,rTrkVtx);
    }
    else if (anaVersion==NuCuts::kCC0325Std ||
             anaVersion == NuCuts::kCC0720Std ||
             anaVersion==NuCuts::kCC0720Test ||
             anaVersion==NuCuts::kJJH1 ||
             anaVersion==NuCuts::kJJE1 ||
             anaVersion==NuCuts::kJJE2 ||
             anaVersion == NuCuts::kNMB0325Alpha ||
             anaVersion == NuCuts::kNMB0325Bravo ||
             anaVersion == NuCuts::kNMB0325BravoTwo ||
             anaVersion == NuCuts::kNMB0325Charlie ||
             anaVersion == NuCuts::kNMB0325Delta ||
             anaVersion == NuCuts::kNMB0325Echo ||
             anaVersion == NuCuts::kRM1 ||
             anaVersion == NuCuts::kRM2 ||
             anaVersion == NuCuts::kNMBFree ||
             anaVersion == NuCuts::kRHC ||
             anaVersion == NuCuts::kMSRock) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on FD xyz point being in kCC0325Std fid vol"<<endl;
      if (!infid_initialised) {
        infid_initialised=true;
        choose_infid_set("cc2008");
      }
      Bool_t isin = infid(static_cast<Detector::Detector_t>(detector),
        static_cast<SimFlag::SimFlag_t>(simFlag),
        x,y,z);
      // Rock muons use an anti-fiducial cut
      if (anaVersion == NuCuts::kMSRock) {
        MAXMSG("NuCuts",Msg::kInfo,1)
          <<"Cutting on INVERSE fiducial volume"<<endl;
        return !isin;
      }
      return isin;
    }
    else if (anaVersion==NuCuts::kFullDST) {
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"For kFullDST cutting on FD xyz point being in"
        <<" \"loose\" fid vol"<<endl;
      return this->IsInFidVolLoose(x,y,z,detector);
    }
    else {
      //use kCC0093Std here since all the variables are safely defined
      MAXMSG("NuCuts",Msg::kInfo,1)
        <<"Cutting on FD xyz point being in kCC0093Std fid vol"<<endl;
      return this->IsInFidVolFDCC0093Std(x,y,z);
    }
  }
  else {
    cout<<"Ahhh, detector not recognised, det="<<detector<<endl;
    return false;
  }
}

//......................................................................

Bool_t NuCuts::IsInFidVolTrueEvt(const NuMCEvent& mc) const
{
  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on IsInFidVolTrueEvt()"<<endl;
  return this->IsInFidVol
    (mc.vtxxMC,mc.vtxyMC,mc.vtxzMC,mc.vtxuMC, mc.vtxvMC,
     mc.planeTrkVtxMC,mc.rTrkVtxMC,
     mc.detector,mc.anaVersion,mc.releaseType,mc.simFlag);
}

//......................................................................

Bool_t NuCuts::IsInFidVolTrueEvt(const NuEvent& nu) const
{
  //MAXMSG("NuCuts",Msg::kInfo,1) <<"Cutting on IsInFidVolTrueEvt()"<<endl;
  if (nu.simFlag != SimFlag::kMC) {
    MSG("NuCuts",Msg::kError) << "Cannot do true cut on data" << endl;
    MSG("NuCuts",Msg::kFatal) << "Cannot do true cut on data" << endl;
  }
  // Not sure if this needs to handle rock muons - but apply anyway
  return this->IsInFidVol
    (nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,nu.vtxuMC, nu.vtxvMC,
     nu.planeTrkVtxMC,nu.rTrkVtxMC,
     nu.detector,nu.anaVersion,nu.releaseType,nu.simFlag);
}

//......................................................................

Bool_t NuCuts::IsInFidVolTrk(const NuEvent& nu) const
{

  //by default use the trk vtx as is
  Float_t z=nu.zTrkVtx;

  //move trk vtx upstream by 3.92cm from scintillator to steel
  //for these analysis versions
  if (nu.anaVersion==NuCuts::kCC0325Std ||
      nu.anaVersion==NuCuts::kCC0720Std ||
      nu.anaVersion==NuCuts::kCC0720Test ||
      nu.anaVersion==NuCuts::kJJH1 ||
      nu.anaVersion==NuCuts::kJJE1 ||
      IsNMB(nu) ) {
    z=nu.zTrkVtx-(0.0392*Munits::m);
  }

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on IsInFidVolTrk()"<<endl;
  return this->IsInFidVol(nu.xTrkVtx,nu.yTrkVtx,
                          z,
                          nu.uTrkVtx,nu.vTrkVtx,
                          nu.planeTrkVtx,nu.rTrkVtx,
                          nu.detector,nu.anaVersion,
                          nu.releaseType,nu.simFlag);
}

//......................................................................

Bool_t NuCuts::IsInFidVolEvt(const NuEvent& nu) const
{

  //by default use the evt vtx as is
  Float_t z=nu.zEvtVtx;

  // TODO CJB I don't think we need to move the vertex, but I'm not sure
  /*
  //move trk vtx upstream by 3.92cm from scintillator to steel
  //for these analysis versions
  if (nu.anaVersion==NuCuts::kCC0325Std ||
      nu.anaVersion==NuCuts::kCC0720Std ||
      nu.anaVersion==NuCuts::kJJH1 ||
      nu.anaVersion==NuCuts::kJJE1 ||
      IsNMB(nu) ) {
    z=nu.zTrkVtx-(0.0392*Munits::m);
  }
  */

  MAXMSG("NuCuts",Msg::kInfo,1)
    <<"Cutting on IsInFidVolTrk()"<<endl;
  return this->IsInFidVol(nu.xEvtVtx,nu.yEvtVtx,
                          z,
                          nu.uEvtVtx,nu.vEvtVtx,
                          nu.planeEvtVtx,nu.rEvtVtx,
                          nu.detector,nu.anaVersion,
                          nu.releaseType,nu.simFlag);
}

//......................................................................

Bool_t NuCuts::IsInFidVolOffset(Float_t x,Float_t y,Float_t z) const
{
  const Float_t zerox=0.9*(Munits::m);
  const Float_t zeroy=0.0*(Munits::m);
  const Float_t minZ=1.0*(Munits::m);
  const Float_t maxZ=4.0*(Munits::m);
  const Float_t maxR=0.6*(Munits::m);
  Bool_t passFid=this->IsInCylindricalVolume(x,y,z,zerox,zeroy,
                                             minZ,maxZ,maxR);
  passFid=passFid && x>=0.7;
  return passFid;
}

//......................................................................

Bool_t NuCuts::IsInFidVolOffset(const NuEvent& nu) const
{
  return this->IsInFidVolOffset(nu.xTrkVtx,nu.yTrkVtx,nu.zTrkVtx);
}

//......................................................................

Bool_t NuCuts::IsInFidVolNDNuMuBar(Float_t x,Float_t y,
                                   Float_t z) const
{
  //fiducial criteria on vtx for near detector
  const Float_t beamzerox=1.4828*(Munits::m);//new for 2.5 analysis
  const Float_t beamzeroy=0.2384*(Munits::m);//new for 2.5 analysis

  //define this volume
  Float_t minZ=0.6*(Munits::m);
  Float_t maxZ=4.0*(Munits::m);
  Float_t maxR=0.8*(Munits::m);

  return this->IsInCylindricalVolume(x,y,z,beamzerox,beamzeroy,
                                     minZ,maxZ,maxR);
}

//......................................................................

Bool_t NuCuts::IsInFidVolNDCC0093Std(Float_t x,Float_t y,
                                     Float_t z) const
{
  //got these numbers from Masaki's spectrum validation talks
  //fiducial criteria on vtx for near detector
  const Float_t beamzerox=1.4885*(Munits::m);
  const Float_t beamzeroy=0.1397*(Munits::m);
  //Float_t beamzeroy = 0.1397-reco_vtxz*TMath::Tan(.058);

  //define this volume
  Float_t minZ=1.0*(Munits::m);
  Float_t maxZ=5.0*(Munits::m);
  Float_t maxR=1.0*(Munits::m);

  return this->IsInCylindricalVolume(x,y,z,beamzerox,beamzeroy,
                                     minZ,maxZ,maxR);
}

//......................................................................

Bool_t NuCuts::IsInFidVolNDCC0250Std(Float_t x,Float_t y,
                                     Float_t z) const
{
  //fiducial criteria on vtx for near detector
  const Float_t beamzerox=1.4828*(Munits::m);//new for 2.5 analysis
  const Float_t beamzeroy=0.2384*(Munits::m);//new for 2.5 analysis

  //define this volume
  Float_t minZ=1.0*(Munits::m);
  Float_t maxZ=5.0*(Munits::m);
  Float_t maxR=1.0*(Munits::m);

  return this->IsInCylindricalVolume(x,y,z,beamzerox,beamzeroy,
                                     minZ,maxZ,maxR);
}

//......................................................................

Bool_t NuCuts::IsInFidVolFDNuMuBar(Int_t planeTrkVtx,
                                   Float_t rTrkVtx) const
{
  //fiducial criteria on vtx for far detector
  //these are the CC-analysis 2.5e20 cuts
  //plus a coil hole cut based on John Marshall's talk DocDB 3644

  //minos-doc-3156
  if (rTrkVtx>(0.4*Munits::m) &&
      rTrkVtx<sqrt(14.)*Munits::m &&
      ((planeTrkVtx>4 && planeTrkVtx<241) ||
       (planeTrkVtx>253 && planeTrkVtx<466))) return true;
  else return false;
}

//......................................................................

Bool_t NuCuts::IsInFidVolFDCC0093Std(Float_t x,Float_t y,
                                     Float_t z) const
{
  //got these numbers from Masaki's spectrum validation talks
  //there is no coil cut!
  //fiducial criteria on vtx for far detector
  if(((z>0.50 && z<14.3) || (z>16.2 && z<28.0)) &&
     ((x*x)+(y*y))<14.0) return true;
  else return false;
}

//......................................................................

Bool_t NuCuts::IsInFidVolFDCC0250Std(Int_t planeTrkVtx,
                                     Float_t rTrkVtx) const
{
  //fiducial criteria on vtx for far detector
  //these are Andy Blake's cuts in DocDB 2728
  //but without the coil cut

  //if (drTrkFidvtx>(0.2*Munits::m) &&
  //    rTrkVtx>(0.4*Munits::m) &&
  //    ((planeTrkVtx>4 && planeTrkVtx<241) ||
  //     (planeTrkVtx>253 && planeTrkVtx<466))) return true;

  //minos-doc-3156
  if (rTrkVtx<sqrt(14.)*Munits::m &&
      ((planeTrkVtx>4 && planeTrkVtx<241) ||
       (planeTrkVtx>253 && planeTrkVtx<466))) return true;
  else return false;
}

//fiducial criteria on vtx for far detector
//these are Andy Blake's cuts in DocDB 2728
//coil cut added wrt 2.5e20 analysis
//and coded using z position instead of plane

//const Float_t minZSM1=0.215*Munits::m;
//const Float_t maxZSM1=14.255*Munits::m;
//const Float_t minZSM2=16.122*Munits::m;
//const Float_t maxZSM2=28.722*Munits::m;
//const Float_t minR=0.4*Munits::m;
//const Float_t maxR=sqrt(14.)*Munits::m;

//calc radius
//Float_t r=sqrt((x*x)+(y*y));

//make the cut
//if ((z>minZSM1 && z<maxZSM1) ||
//      (z>minZSM2 && z<maxZSM2) &&
//      r>minR && r<maxR) return true;
//else return false;

//......................................................................

Bool_t NuCuts::IsInFidVolPitt(Float_t x,Float_t y,Float_t z,
                              Float_t u,Float_t v) const
{
  // Is the event vertex, (x,y,u,v,z) inside the fiducial volume
  // used by the pittsburgh group?
  // code from Mike
  if( !( z>0.6 && z<3.56) ) return kFALSE;
  if( !( u>0.3 && u<1.8) ) return kFALSE;
  if( !( v>-1.8 && v< -0.3) ) return kFALSE;
  if( x>=2.4) return kFALSE;
  const Float_t coil_cut=0.8*0.8;
  if( (x*x + y*y) <= coil_cut) return kFALSE;
  return kTRUE;
}

//......................................................................

Bool_t NuCuts::IsInFidVolLoose(const NuEvent& nu) const
{

  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  Bool_t trkInFidLoose1=false;
  if (nu.trkExists1) trkInFidLoose1=lib.cuts.IsInFidVolLoose
                       (nu.xTrkVtx1,nu.yTrkVtx1,nu.zTrkVtx1,
                        nu.detector);
  Bool_t trkInFidLoose2=false;
  if (nu.trkExists2) trkInFidLoose2=lib.cuts.IsInFidVolLoose
                       (nu.xTrkVtx2,nu.yTrkVtx2,nu.zTrkVtx2,
                        nu.detector);
  Bool_t trkInFidLoose3=false;
  if (nu.trkExists3) trkInFidLoose3=lib.cuts.IsInFidVolLoose
                       (nu.xTrkVtx3,nu.yTrkVtx3,nu.zTrkVtx3,
                        nu.detector);
  Bool_t shwInFidLoose1=false;
  if (nu.shwExists1) shwInFidLoose1=lib.cuts.IsInFidVolLoose
                       (nu.xShwVtx1,nu.yShwVtx1,nu.zShwVtx1,
                        nu.detector);
  Bool_t shwInFidLoose2=false;
  if (nu.shwExists2) shwInFidLoose2=lib.cuts.IsInFidVolLoose
                       (nu.xShwVtx2,nu.yShwVtx2,nu.zShwVtx2,
                        nu.detector);
  Bool_t shwInFidLoose3=false;
  if (nu.shwExists3) shwInFidLoose3=lib.cuts.IsInFidVolLoose
                       (nu.xShwVtx3,nu.yShwVtx3,nu.zShwVtx3,
                        nu.detector);
  const Bool_t trkInFidLoose=trkInFidLoose1 || trkInFidLoose2 ||
    trkInFidLoose3;
  const Bool_t shwInFidLoose=shwInFidLoose1 || shwInFidLoose2 ||
    shwInFidLoose3;
  const Bool_t evtInFidLoose=lib.cuts.IsInFidVolLoose
    (nu.xEvtVtx,nu.yEvtVtx,nu.zEvtVtx,nu.detector);

  //cut out events that don't have either trk or shw or evt
  //in loose fid vol
  return (trkInFidLoose || shwInFidLoose || evtInFidLoose);
}

//......................................................................

Bool_t NuCuts::IsInFidVolLoose(Float_t x,Float_t y,Float_t z,
                               Int_t detector) const
{

  //return true for all FD events since almost all are in FV anyway
  if (detector==Detector::kFar) return true;

  //use the 0093CC analysis beam centre since it is more central
  //than the 0250CC analysis
  const Float_t beamzerox=1.4828*(Munits::m);
  const Float_t beamzeroy=0.1397*(Munits::m);
  const Float_t minZ=0.5*(Munits::m);
  const Float_t maxZ=5.1*(Munits::m);
  const Float_t maxR=1.15*(Munits::m);//allows 1m for zero-y & bc=0.23

  //check if in volume
  Bool_t passFid=this->IsInCylindricalVolume
    (x,y,z,beamzerox,beamzeroy,minZ,maxZ,maxR);

  //return selection
  return passFid;
}

//......................................................................

Bool_t NuCuts::IsInCylindricalVolume(Float_t x,Float_t y,
                                     Float_t z,
                                     Float_t zerox,Float_t zeroy,
                                     Float_t minZ,Float_t maxZ,
                                     Float_t maxR) const
{
  Float_t r=sqrt(pow((x-zerox),2)+pow((y-zeroy),2));
  if(z>minZ && z<maxZ && r<maxR) return true;
  else return false;
}

//......................................................................

void NuCuts::FindZCuts() const
{

  //NOTE: using trk.vtx.z-0.0392 for CC analysis to move vtx from scint
  //to steel where interaction happened

  //steel planes 14-68 inclusive for CC2008 analysis
  //ND LOWER CUT: actually want cut to be in the large air gap
  //upstream of plane 14, so use 13
  //ND UPPER CUT: want cut to be in large air gap downstream
  //of 68, so use 68

  //FD SM1 LOWER CUT: use plane 3 (based on trk.vtx.z-0.0392)
  //FD SM2 UPPER CUT: use plane 239
  //FD SM1 LOWER CUT: use plane 252
  //FD SM2 UPPER CUT: use plane 464

  //ND: trk.vtx.z-0.0392>0.809 && trk.vtx.z-0.0392<4.076

  //(trk.vtx.z-0.0392>0.215 && trk.vtx.z-0.0392<14.255) ||                      (trk.vtx.z-0.0392>16.122 && trk.vtx.z-0.0392<28.722)
  //These positions are based on MC (FD data is shifted by ~0.7-1.4cm wrt to MC)
  //Andy's cuts are currently in plane form:
  //(trk.vtx.plane>4 && trk.vtx.plane<241) || (trk.vtx.plane>253 && trk.vtx.plane<466)

  //enum EAirRegion {
  //kBetweenScintAndSteel =  0,
  //kDownstreamOfSteel          = +1
  //};
  //enum EPosition {
  //kUpstream           = -1,
  //kMiddle                     =  0,
  //kDownstream = +1
  //};

  FidVol::find_z_cuts(17,84,8,240,255,452,
                      FidVol::kBetweenScintAndSteel,
                      FidVol::kMiddle,Munits::micrometer);

  FidVol::find_z_cuts(17,84,8,240,255,452,
                      FidVol::kDownstreamOfSteel,
                      FidVol::kMiddle,Munits::micrometer);

  //use CC2008 cuts
  FidVol::find_z_cuts(13,68,3,239,252,464,
                      FidVol::kDownstreamOfSteel,
                      FidVol::kMiddle,Munits::micrometer);
}

//......................................................................

Bool_t NuCuts::IsGoodUVVtx(const NuEvent& nu) const
{
  Float_t cutValue=5;
  if (false) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on UV trk vtx<="<<cutValue<<endl;
    if (abs(nu.trkVtxUVDiffPl)<=cutValue) return true;
    else return false;
  }
  else if (nu.anaVersion==NuCuts::kCC0093Std ||
           nu.anaVersion==NuCuts::kCC0250Std ||
           nu.anaVersion==NuCuts::kJJE1 ||
           nu.anaVersion==NuCuts::kCC0325Std ||
           nu.anaVersion==NuCuts::kCC0720Std ||
           nu.anaVersion==NuCuts::kCC0720Test ||
           nu.anaVersion==NuCuts::kJJH1 || IsNMB(nu) ) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Not cutting on UV trk vtx"<<endl;
    return true;
  }
  else {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Cutting on UV trk vtx<="<<cutValue<<endl;
    if (abs(nu.trkVtxUVDiffPl)<=cutValue) return true;
    else return false;
  }
}

//......................................................................
/*
  void NuCuts::FillTimeHistos(const NtpStRecord& ntp,
  const NtpSREvent& evt) const
  {

  static TH1F* hTimeLength=0;
  static TH1F* hTimeBefore=0;
  static TH1F* hTimeAfter=0;

  static TH1F* hSlcTimeLength=0;
  static TH1F* hSlcTimeBefore=0;
  static TH1F* hSlcTimeAfter=0;

  static TH1F* hTrkTimeLength=0;
  static TH1F* hTrkTimeBefore=0;
  static TH1F* hTrkTimeAfter=0;

  static TH1F* hMedianTrkSlcDiff=0;

  if (hTimeLength==0) {
  MSG("NuCuts",Msg::kInfo)
  <<"Creating time histos..."<<endl;

  hTimeLength=new TH1F("hTimeLength","hTimeLength",10000,-100,50000);
  hTimeLength->SetFillColor(0);
  hTimeLength->SetTitle("Snarl Length in Time");
  hTimeLength->GetXaxis()->SetTitle("Time (ns)");
  hTimeLength->GetXaxis()->CenterTitle();

  hTimeBefore=new TH1F("hTimeBefore","hTimeBefore",10000,-100,50000);
  hTimeBefore->SetFillColor(0);
  hTimeBefore->SetTitle("Time Between Event and First Hit");
  hTimeBefore->GetXaxis()->SetTitle("Time (ns)");
  hTimeBefore->GetXaxis()->CenterTitle();

  hTimeAfter=new TH1F("hTimeAfter","hTimeAfter",10000,-100,50000);
  hTimeAfter->SetFillColor(0);
  hTimeAfter->SetTitle("Time Between Event and Last Hit");
  hTimeAfter->GetXaxis()->SetTitle("Time (ns)");
  hTimeAfter->GetXaxis()->CenterTitle();

  hSlcTimeLength=new TH1F
  ("hSlcTimeLength","hSlcTimeLength",10000,-100,50000);
  hSlcTimeLength->SetFillColor(0);
  hSlcTimeLength->SetTitle("Slice Length in Time");
  hSlcTimeLength->GetXaxis()->SetTitle("Time (ns)");
  hSlcTimeLength->GetXaxis()->CenterTitle();

  hSlcTimeBefore=new TH1F
  ("hSlcTimeBefore","hSlcTimeBefore",10000,-100,50000);
  hSlcTimeBefore->SetFillColor(0);
  hSlcTimeBefore->SetTitle("Time Between Event and First Hit in Slice");
  hSlcTimeBefore->GetXaxis()->SetTitle("Time (ns)");
  hSlcTimeBefore->GetXaxis()->CenterTitle();

  hSlcTimeAfter=new TH1F
  ("hSlcTimeAfter","hSlcTimeAfter",10000,-100,50000);
  hSlcTimeAfter->SetFillColor(0);
  hSlcTimeAfter->SetTitle("Time Between Event and Last Hit in Slice");
  hSlcTimeAfter->GetXaxis()->SetTitle("Time (ns)");
  hSlcTimeAfter->GetXaxis()->CenterTitle();

  hTrkTimeLength=new TH1F
  ("hTrkTimeLength","hTrkTimeLength",10000,-100,50000);
  hTrkTimeLength->SetFillColor(0);
  hTrkTimeLength->SetTitle("Track Length in Time");
  hTrkTimeLength->GetXaxis()->SetTitle("Time (ns)");
  hTrkTimeLength->GetXaxis()->CenterTitle();

  hTrkTimeBefore=new TH1F
  ("hTrkTimeBefore","hTrkTimeBefore",10000,-100,50000);
  hTrkTimeBefore->SetFillColor(0);
  hTrkTimeBefore->SetTitle("Time Between Event and First Tracked Hit");
  hTrkTimeBefore->GetXaxis()->SetTitle("Time (ns)");
  hTrkTimeBefore->GetXaxis()->CenterTitle();

  hTrkTimeAfter=new TH1F
  ("hTrkTimeAfter","hTrkTimeAfter",10000,-100,50000);
  hTrkTimeAfter->SetFillColor(0);
  hTrkTimeAfter->SetTitle("Time Between Event and Last Tracked Hit");
  hTrkTimeAfter->GetXaxis()->SetTitle("Time (ns)");
  hTrkTimeAfter->GetXaxis()->CenterTitle();

  hMedianTrkSlcDiff=new TH1F
  ("hMedianTrkSlcDiff","hMedianTrkSlcDiff",10000,-3000,3000);
  hMedianTrkSlcDiff->SetFillColor(0);
  hMedianTrkSlcDiff->SetTitle("Diff. between Median Slc and Median Trk Time");
  hMedianTrkSlcDiff->GetXaxis()->SetTitle("Time (ns)");
  hMedianTrkSlcDiff->GetXaxis()->CenterTitle();
  }

  TClonesArray& stpTca=(*ntp.stp);
  Int_t numStps=stpTca.GetEntries();

  Double_t earliestTime=999;
  Double_t latestTime=-999;
  multiset<Double_t> times;

  //this loop is over ALL the strips
  for (Int_t i=0;i<numStps;i++){

  const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(stpTca[i]);
  const NtpSRStrip& stp=(*pstp);

  //Int_t pl=stp.plane;

  //get earliest and latest times in snarl
  if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
  if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
  if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=stp.time1;
  if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=stp.time0;

  Double_t time=stp.time0;
  if (time<0 || time>1) time=stp.time1;

  if (time>0 && time<=1) {
  times.insert(time);
  }

  MAXMSG("NuCuts",Msg::kDebug,500)
  <<"t0="<<stp.time0<<", t1="<<stp.time1
  <<", lT="<<latestTime<<", eT="<<earliestTime
  <<", dT="<<(latestTime-earliestTime)/Munits::ns<<endl;
  }

  //get the median time from the map
  multiset<Double_t>::const_iterator it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTime=*it;
  MAXMSG("NuCuts",Msg::kDebug,200)
  <<"Median time="<<medianTime
  <<", dBefore="<<(medianTime-earliestTime)/Munits::ns
  <<", dAfter="<<(latestTime-medianTime)/Munits::ns<<endl;

  hTimeAfter->Fill((latestTime-medianTime)/Munits::ns);
  hTimeBefore->Fill((medianTime-earliestTime)/Munits::ns);
  hTimeLength->Fill((latestTime-earliestTime)/Munits::ns);

  //now look at slice times

  //reset variables
  earliestTime=999;
  latestTime=-999;
  times.clear();

  TClonesArray& slcTca=(*ntp.slc);

  //get the slice associated with this event
  const NtpSRSlice* pslc=
  dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
  const NtpSRSlice& slc=(*pslc);

  //loop over strips in slc
  for (Int_t i=0;i<slc.nstrip;++i) {
  const NtpSRStrip* pstp=
  dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
  const NtpSRStrip& stp=(*pstp);

  //get earliest and latest times in snarl
  if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
  if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
  if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=stp.time1;
  if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=stp.time0;

  Double_t time=stp.time0;
  if (time<0 || time>1) time=stp.time1;

  if (time>0 && time<=1) {
  times.insert(time);
  }
  }

  //get the median time from the map
  it=times.begin();
  advance(it,times.size()/2);
  medianTime=*it;

  hSlcTimeAfter->Fill((latestTime-medianTime)/Munits::ns);
  hSlcTimeBefore->Fill((medianTime-earliestTime)/Munits::ns);
  hSlcTimeLength->Fill((latestTime-earliestTime)/Munits::ns);

  //now look at track times

  //reset variables
  earliestTime=999;
  latestTime=-999;
  times.clear();

  TClonesArray& trkTca=(*ntp.trk);
  Int_t numTrks=evt.ntrack;

  for (Int_t itrk=0;itrk<numTrks;itrk++){
  const NtpSRTrack* ptrk=
  dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
  const NtpSRTrack& trk=(*ptrk);

  TClonesArray& stpTca=(*ntp.stp);

  //this loop is just over the tracked strips
  for (Int_t i=0;i<trk.nstrip;i++){

  const NtpSRStrip* pstp=
  dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
  const NtpSRStrip& stp=(*pstp);

  //get earliest and latest times in snarl
  if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
  if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
  if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=
  stp.time1;
  if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=
  stp.time0;

  Double_t time=stp.time0;
  if (time<0 || time>1) time=stp.time1;

  if (time>0 && time<=1) {
  times.insert(time);
  }
  }
  }

  //get the median time from the map
  it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTimeTrk=*it;

  MAXMSG("NuCuts",Msg::kDebug,500)
  <<"Median time for slc and trk difference="
  <<(medianTime-medianTimeTrk)/Munits::ns<<endl;

  hMedianTrkSlcDiff->Fill((medianTime-medianTimeTrk)/Munits::ns);

  hTrkTimeAfter->Fill((latestTime-medianTimeTrk)/Munits::ns);
  hTrkTimeBefore->Fill((medianTimeTrk-earliestTime)/Munits::ns);
  hTrkTimeLength->Fill((latestTime-earliestTime)/Munits::ns);
  }
*/
//......................................................................
/*
  void NuCuts::PrintNtpSt(const NtpStRecord& ntp) const
  {
  TClonesArray& tcaTk=(*ntp.trk);
  Int_t numTrks=tcaTk.GetEntries();
  TClonesArray& stpTca=(*ntp.stp);
  Int_t numStps=stpTca.GetEntries();

  MAXMSG("NuCuts",Msg::kInfo,200)
  <<"PrintNtpSt: numTrks="<<numTrks
  <<", total numStps="<<numStps<<endl;

  //Loop over tracks
  for (Int_t itrk=0;itrk<numTrks;itrk++){
  const NtpSRTrack* ptrk=
  dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
  const NtpSRTrack& trk=(*ptrk);

  TClonesArray& stpTca=(*ntp.stp);

  //this loop is just over the tracked strips
  for (Int_t i=0;i<trk.nstrip;i++){

  const NtpSRStrip* pstp=
  dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
  const NtpSRStrip& stp=(*pstp);
  MAXMSG("NDNuCuts",Msg::kInfo,20000)
  <<"Trk: st="<<stp.strip<<", pl="<<stp.plane
  <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
  }
  }

  //this loop is over ALL the strips
  for (Int_t i=0;i<numStps;i++){
  const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(stpTca[i]);
  const NtpSRStrip& stp=(*pstp);

  MAXMSG("NDNuCuts",Msg::kInfo,20000)
  <<"All strips: st="<<stp.strip<<", pl="<<stp.plane
  <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
  }
  }
*/
//......................................................................

void NuCuts::CalcTotalPot(Float_t& totalPotHist,
                          Float_t& totalPotBadHist) const
{
  //calc the total POT
  TH1F* hPottortgt=dynamic_cast<TH1F*>(gROOT->FindObject("hPottortgt"));
  TH1F* hPotBadtortgt=
    dynamic_cast<TH1F*>(gROOT->FindObject("hPotBadtortgt"));

  if (hPottortgt && hPotBadtortgt) {
    //hPottortgt->Print();
    //hPotBadtortgt->Print();

    for (Int_t i=0;i<hPotBadtortgt->GetNbinsX();i++) {
      Float_t num=hPottortgt->GetBinContent(i);
      Float_t val=hPottortgt->GetBinLowEdge(i);
      val+=hPottortgt->GetBinWidth(i)/2.;
      Float_t sum=num*val;
      if (sum<0) {
        MAXMSG("NuCuts",Msg::kInfo,10)
          <<"sum<0, sum="<<sum
          <<", num="<<num
          <<", val="<<val<<endl;
        sum=0;
      }
      totalPotHist+=sum;

      Float_t numBad=hPotBadtortgt->GetBinContent(i);
      Float_t valBad=hPotBadtortgt->GetBinLowEdge(i);
      valBad+=hPotBadtortgt->GetBinWidth(i)/2.;
      Float_t sumBad=numBad*valBad;
      if (sumBad<0) {
        MAXMSG("NuCuts",Msg::kInfo,10)
          <<"sumBad<0, sumBad="<<sumBad
          <<", num="<<numBad
          <<", val="<<valBad<<endl;
        sumBad=0;
      }
      totalPotBadHist+=sumBad;

      MSG("NuCuts",Msg::kDebug)
        <<"i="<<i<<", Good: num="<<num
        <<", val="<<val<<", sum="<<sum
        <<", totalPot="<<totalPotHist<<endl;

      MSG("NuCuts",Msg::kDebug)
        <<"i="<<i<<", Bad:      num="<<numBad
        <<", val="<<valBad<<", sum="<<sumBad
        <<", totalPotBad="<<totalPotBadHist<<endl;
    }
  }
  else cout<<"Can't find histograms to calc total POT"<<endl;
}

//......................................................................

void NuCuts::CheckTrackDirectionIsPositive(const NuEvent& nu) const
{
  if (nu.planeTrkEnd-nu.planeTrkVtx<0) {
    //NuPlots plots;
    //plots.PrintRunSnarlEvent(cout,nu);
    MSG("NuAnalysis",Msg::kWarning)
      <<endl<<endl<<endl
      <<"Track direction wrong!!!, vtx plane="<<nu.planeTrkVtx
      <<", end plane="<<nu.planeTrkEnd
      <<endl<<endl<<endl<<endl;
  }
}

//......................................................................

const Char_t* NuCuts::AsString(NuAnaVersion_t v) const
{
  // return char string
  switch (v) {
  case kCC0093Std: return "Std CC 1.27e20 POT";
  case kJJH1:      return "JJH NuMuBar";
  case kJJE1:      return "JJE NuMuBar";
  case kJJE2:      return "JJE thesis NuMuBar";
  case kCC0250Std: return "Std CC 2.50e20 POT";
  case kCC0325Std: return "Std CC 3.25e20 POT";
  case kCC0720Test: return "Std CC (with JMCUTS) for 2010 analysis";
  case kFullDST:   return "Full DST w/ loose cuts";
  case kMSRock:    return "Matt Straits Ely2009 Rock muon selection";
  case kCC0720Std: return "Std CC 7.20e20 POT";
  case kUnknown:   return "Unknown"; break;
  default:         return "!?Bad Value Passed?!"; break;
  }
}

//......................................................................

const Double_t NuCuts::FiducialMass(const NuEvent& nu) const
{
  return this->FiducialMass(nu.detector,nu.simFlag,nu.anaVersion);
}

//......................................................................

const Double_t NuCuts::FiducialMass(const Int_t detector,
                                    const Int_t simFlag,
                                    const Int_t anaVersion) const
{
  if (detector == Detector::kNear){
    if (anaVersion==NuCuts::kJJH1 ||
        anaVersion==NuCuts::kJJE1 ||
        anaVersion==NuCuts::kJJE2 ||
        anaVersion == NuCuts::kCC0325Std ||
        anaVersion == NuCuts::kCC0720Std ||
        anaVersion == NuCuts::kNMB0325Alpha ||
        anaVersion == NuCuts::kNMB0325Bravo ||
        anaVersion == NuCuts::kNMB0325BravoTwo ||
        anaVersion == NuCuts::kNMB0325Charlie ||
        anaVersion == NuCuts::kNMB0325Delta ||
        anaVersion == NuCuts::kNMB0325Echo ||
        anaVersion == NuCuts::kRM1 ||
        anaVersion == NuCuts::kRM2 ||
        anaVersion == NuCuts::kNMBFree ||
        anaVersion == NuCuts::kRHC){
      //return 23.4200178*Munits::tonne;//JJE
      if (simFlag==SimFlag::kData) return 23.7217*
                                     Munits::tonne;//JJH
      else if (simFlag==SimFlag::kMC) return 23.4989*
                                        Munits::tonne;//JJH
      else {
        cout<<"Ahhhh simFlag="<<simFlag<<endl;
        return -1;
      }
    }
    else if (anaVersion == NuCuts::kCC0093Std){
      return -1.0;
    }
    else if (anaVersion == NuCuts::kCC0250Std){
      return 44.66068761*Munits::tonne;
    }
    else {
      //PittFidVol
      return -1.0;
    }
  }
  else if (detector == Detector::kFar){
    if (anaVersion==NuCuts::kJJH1 ||
        anaVersion==NuCuts::kJJE1 ||
        anaVersion==NuCuts::kJJE2 ||
        anaVersion == NuCuts::kCC0325Std ||
        anaVersion == NuCuts::kCC0720Std ||
        anaVersion == NuCuts::kCC0720Test ||
        anaVersion == NuCuts::kNMB0325Alpha ||
        anaVersion == NuCuts::kNMB0325Bravo ||
  anaVersion == NuCuts::kNMB0325BravoTwo ||
        anaVersion == NuCuts::kNMB0325Charlie ||
        anaVersion == NuCuts::kNMB0325Delta ||
        anaVersion == NuCuts::kNMB0325Echo ||
        anaVersion == NuCuts::kRM1 ||
        anaVersion == NuCuts::kRM2 ||
        anaVersion == NuCuts::kNMBFree ||
        anaVersion == NuCuts::kRHC) {
      //return 4.119215650*Munits::kilotonne;//JJE
      if (simFlag==SimFlag::kData) return 4.17235*
                                     Munits::kilotonne;//JJH
      else if (simFlag==SimFlag::kMC) return 4.13923*
                                        Munits::kilotonne;//JJH
      else {
        cout<<"Ahhhh simFlag="<<simFlag<<endl;
        return -1;
      }
    }
    else if (anaVersion == NuCuts::kCC0093Std){
      return -1.0;
    }
    else if (anaVersion == NuCuts::kCC0250Std){
      return 4.164776186*Munits::kilotonne;
      //Neglecting the existance of the coil hole!
    }
    else {
      //IsInFidVolFDCC0093Std
      return -1.0;
    }
  }
  else{
    MAXMSG("NuCuts",Msg::kWarning,5)
      << "Bad detector: returning 0 fiducial mass." << endl;
    return 0;
  }
}

//......................................................................
Bool_t NuCuts::FreeCuts(const NuEvent& nu, const NuXMLConfig* xmlConfig)const{
  if(nu.roID<=xmlConfig->RoID()||
     nu.jmID<=xmlConfig->JmID()||
     nu.jmEventknnID <= xmlConfig->NtID()||
     nu.dpID<=xmlConfig->DpID() ||
     nu.poIDKin<=xmlConfig->PoKIN() ||
     fabs(1./nu.sigqp_qp)<=xmlConfig->Sigqp() ||
     fabs(nu.relativeAngle - TMath::Pi())<= xmlConfig->RelativeAngle() ||
     nu.trkLength<=xmlConfig->TrackLength() ||
     nu.prob<=xmlConfig->Trkpro() ||
     nu.smoothMajC<=xmlConfig->MajorityCurvature())
    {
      MSG("NuCuts",Msg::kInfo)
        <<"Inside free cuts true"<<endl;
      return false;}
  else
    {MSG("NuCuts",Msg::kInfo)
      <<"Inside free cuts false"<<endl;
      return true;}
}
//......................................................................

// Rock muon cut variables
Bool_t NuCuts::IsGoodCosPr(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kMSRock) {
    // Rock-like if cosPrTrkVtx < -0.2
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making Cut: IsGoodCosPr < -0.2" << endl;
    if (nu.cosPrTrkVtx >= -0.2) return false;
  }

  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodTrackDr(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kMSRock) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making Cut: IsGoodTrackDr < 0.06" << endl;

    // Rock-like if drTrkFidVtx < 0.06
    if (nu.drTrkFidvtx >= 0.06) return false;
  }

  return true;
}

//......................................................................

Bool_t NuCuts::IsGoodNearShowerEn(const NuEvent& nu) const
{
  if (nu.anaVersion==NuCuts::kMSRock) {
    MAXMSG("NuCuts",Msg::kInfo,1)
      <<"Making Cut: IsGoodNearshowerEnergy < 0.5" << endl;

    // Rock-like if trkShwEnNear < 0.5
    if (nu.trkShwEnNear >= 0.5) return false;
  }

  return true;
}

//......................................................................

Bool_t NuCuts::GoodTimeToNearestSpill(const NuEvent &nu, Double_t lower_us, Double_t upper_us)
{
  // Convert the timings to microseconds
  lower_us *= (Munits::microsecond);
  upper_us *= (Munits::microsecond);

  if (nu.timeEvtMin-nu.timeToNearestSpill>lower_us &&
      nu.timeEvtMin-nu.timeToNearestSpill<upper_us) {
    return true;
  }
  else return false;

}

---