NuCutImps.cxx

Go to the documentation of this file.

// NuCutImplementations
// This holds the actual cut implementations i.e. the actual cuts

#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuCuts.h"

#include "MessageService/MsgService.h"

#include "NuCutImps.h"

#include "Conventions/Munits.h"
#include "Conventions/Detector.h"
#include "Conventions/BeamType.h"
#include "Conventions/SimFlag.h"
#include "DataUtil/infid.h"

#include "TFile.h"
#include "TDirectory.h"
#include "TTree.h"
#include "TMath.h"

CVSID("$Id: NuCutImps.cxx,v 1.24 2010/02/08 18:08:54 bckhouse Exp $");

// Do everything in an embedded namespace, because we don't want to
// pollute the global namespace with cut names
namespace NuCutImps {

// CC 'Standard' selection

  CC0325Std::CC0325Std(const NuPlots *plots) : 
    NuCut("CC0325Std", plots)
  {
    SetFidVol("cc2008");
    SetAnaVersion(6);
  }
  
  // The CC selection uses a slightly altered definition of fiducial volume
  // origin: NuCuts::IsInFidVolTrk circa Sept '09
  Bool_t CC0325Std::InFidVol(const NuEvent &nu) const
  {
    // Move trk vtx upstream by 3.92cm from scintillator to steel
    // for this analysis version
    NuEvent nuc = nu;
    nuc.zTrkVtx = nu.zTrkVtx - (0.0392*Munits::m);
    
    return NuCut::InFidVol(nuc);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void CC0325Std::Preselection(const NuEvent &nu)
  {
    if (nu.simFlag==SimFlag::kData) {
      MAXMSG("CC0325Std",Msg::kWarning, 10)
        << "NuCut derivative CC0325Std has not been validated against data. "
        << "Use at your own risk! (Or validate it then remove this warning!)" << endl;
    }
    // Good tracks
    Keep_If(nu.ntrk >= 1, "GoodNumberofTracks");
    
    // Filter out LI events
    Cut_If(nu.isLI || nu.litime != -1, "LICut");
    
    // Is the coil okay? Only cut if we have the cutOnDataQuality flag
    // Keep_Data_If(fNuCuts.IsCoilOk(nu), nu);
    Keep_Data_If(!nu.cutOnDataQuality || nu.coilIsOk, nu, "IsCoilOkay");
  
    // IsGoodBeamDetPOTCountingStage
    // This can mess up POT counting apparently, so be careful!
    Keep_Data_If(nu.goodBeamSntp && nu.coilIsOk, nu,"GoodBeamDetPOT");
    
    // Cut on good beam (IsGoodBeam)
    Keep_Data_If(!nu.cutOnBeamInfo || IsGoodBeam(nu), nu, "GoodBeam");
    
    // Is it in the fiducial volume?
    Keep_If(InFidVol(nu), "FidVol");
    
    // Per-detector cuts
    if (nu.detector == Detector::kFar)
    {
      // TrackFitPass
      Keep_If(nu.trkfitpass == 1, "trackfitpass");
      
      // Data Quality
      Keep_Data_If(nu.isGoodDataQuality || !nu.cutOnDataQuality, nu, "DataQuality");
      
      // IsGoodTimeToNearestSpill (FD only)
      // CC0325Std Uses -20us, +30us
      Keep_Data_If(NuCuts::GoodTimeToNearestSpill(nu, -20, +30) || !nu.cutOnSpillTiming, nu, "SpillTime");
      
      // Cut on track direction angle
      Keep_If(nu.dirCosNu > 0.6, "dirCosNu"); 
    }
    else if (nu.detector == Detector::kNear)
    {
      // Need to implement ND Trackfitpass (uses reclamation)
      // See Zeynep talk DocDB-6382 - mentions 'Nikis Reclamation'
      Keep_If(NuCuts::IsGoodTrackFitPassReclamation(nu), "trackfitreclamation");
      
      // Cut on the coil current
      Keep_Data_If(nu.coilCurrent < -1000 || !nu.cutOnDataQuality, nu, "CoilCurrent");
      
    } // End of per-detector cuts
  } // Eend of preselection function
  
  // The Main Selection ////////////////////////////////////////////////
  void CC0325Std::Selection(const NuEvent &nu)
  {
    Float_t cutValue=0.3;
    // Cut on RoID
    Keep_If(nu.roID > cutValue, "RoID");
    
  }
  
  Bool_t CC0325Std::IsGoodBeam(const NuEvent &nu)
  {
    // Don't measure this unless we have data
    if (nu.simFlag!=SimFlag::kData) {
      // MAXMSG("NuCuts",Msg::kInfo,1)
      //   <<"Not cutting on GoodBeam for simFlag="<<nu.simFlag<<endl;
      return true;
    }

    // Don't cut on this if we've set not to
    if (!nu.cutOnBeamInfo) {
      MAXMSG("CC0325Std",Msg::kWarning,1)
        <<"Not cutting on BeamInfo (flag set to not cut)"<<endl;
      return true;
    }
    
    // Assume we are good, and start from there
    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("CC0325Std",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("CC0325Std",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("CC0325Std",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("CC0325Std",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("CC0325Std",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..."<<endl;
        }
        return false;
      }
    }
    //Beam consistent with first event in DST, so we can let it pass
    return true;
    
  }
  // 
  
  CC0720Test::CC0720Test(const NuPlots *plots):
    NuCut("CC0325Test",plots),
    fNuCutsSelection(NuCuts::kCC0720Test, plots)
  {
    SetFidVol("cc2008");
  }
                
  void CC0720Test::Preselection(const NuEvent &nu)
  {
    // Use old preselection 
    Defer_Preselection(fNuCutsSelection, nu);

  } // end of preselection function   
  void CC0720Test::Selection(const NuEvent &nu)
  {
    Float_t cutValue1=0.3;
    Float_t cutValue2=0.5;
    Float_t cutValue3=0.8;
    // Cut on RoID                                                                                                   
    Keep_If((nu.roID > cutValue1)||
            (nu.jmID > cutValue2)||
            (nu.jmEventknnID > cutValue3)  , "RoID");

  }
  
// NuBar Bravo Selection
  
  Bravo::Bravo(const NuPlots *plots) :
    NuCut("Bravo", plots), 
    fNuCutsSelection(NuCuts::kNMB0325Bravo, plots)
  {
    SetFidVol("cc2008");
    SetAnaVersion(9);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void Bravo::Preselection(const NuEvent &nu)
  {
    // Use old preselection
    Defer_Preselection(fNuCutsSelection, nu);

  } // Eend of preselection function
  
  // The Main Selection ////////////////////////////////////////////////
  void Bravo::Selection(const NuEvent &nu)
  {
    if (nu.charge != +1) { 
      Defer_Selection(fCC0325Std, nu);
    }
    else {
      Keep_If(nu.dpID > 0.25, "DpID");
      Keep_If(1./nu.sigqp_qp > 3.5, "SigmaQP_QP");
      Keep_If(TMath::Abs(nu.relativeAngle - TMath::Pi()) > 2.12, "RelativeAngle");
    }
  }

// NuBar Bravo Selection
  
  BravoPrime::BravoPrime(const NuPlots *plots) : 
    NuCut("BravoPrime", plots), 
    fNuCutsSelection(NuCuts::kNMB0325Bravo, plots)
  {
    SetFidVol("cc2008");
    SetAnaVersion(9);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void BravoPrime::Preselection(const NuEvent &nu)
  {
    // Use old preselection
    Defer_Preselection(fNuCutsSelection, nu);

  } // Eend of preselection function
  
  // The Main Selection ////////////////////////////////////////////////
  void BravoPrime::Selection(const NuEvent &nu)
  {
    if (nu.charge != +1) { 
      Defer_Selection(fCC0325Std, nu);
    }
    else {
      Keep_If(nu.dpID > 0.25, "DpID");
      Keep_If(1./nu.sigqp_qp > 3.6, "SigmaQP_QP");
      Keep_If(TMath::Abs(nu.relativeAngle - TMath::Pi()) > 2.14, "RelativeAngle");
    }
  }
  
  
  // NuBar Charlie Selection
  
  Charlie::Charlie(const NuPlots *plots) :
    NuCut("Charlie", plots),
    fNuCutsSelection(NuCuts::kNMB0325Charlie, plots)
  {
    SetFidVol("cc2008");
    SetAnaVersion(8);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void Charlie::Preselection(const NuEvent &nu)
  {
    // Use old preselection
    Defer_Preselection(fNuCutsSelection, nu);

  }
  
  // The Main Selection ////////////////////////////////////////////////
  void Charlie::Selection(const NuEvent &nu)
  {
    Float_t cutValue=0.625;
    if (nu.charge != +1) { 
      Defer_Selection(fCC0325Std, nu);
    }
    else {
      Keep_If(nu.roID > cutValue, "RoID");
      Keep_If(nu.smoothMajC > -0.22, "MajorityCurv");
      Keep_If(nu.trkLength > 35, "TrackLength");
    }
  }
  
  
  
  
// NuBar Delta Selection
  
  Delta::Delta(const NuPlots *plots) : 
    NuCut("Delta", plots),
    fNuCutsSelection(NuCuts::kNMB0325Delta, plots)
  {
    SetFidVol("cc2008");
    SetAnaVersion(7);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void Delta::Preselection(const NuEvent &nu)
  {
    // Use old preselection
    Defer_Preselection(fNuCutsSelection, nu);

  } // Eend of preselection function
  
  // The Main Selection ////////////////////////////////////////////////
  void Delta::Selection(const NuEvent &nu)
  {
    Float_t cutValue=0.55;
    if (nu.charge != +1) { 
      Defer_Selection(fCC0325Std, nu);
    }
    else {
      Keep_If(nu.roID > cutValue, "RoID");
      Keep_If(nu.prob > 0.003, "FitProb");
      Keep_If(nu.trkLength > 37, "TrackLength");
    }
  }
  
// A common Preselection for CC and NC in the CC2010 analysis

  CCAPresel::CCAPresel(const NuPlots *plots) : 
    NuCut("CCAPresel", plots)
  {

  }
  
  // The Preselection //////////////////////////////////////////////////
  void CCAPresel::Preselection(const NuEvent &nu)
  {
    // IsGoodBeamDetPOTCountingStage
    // This can mess up POT counting apparently, so be careful!
    Keep_Data_If(nu.goodBeamSntp && nu.coilIsOk, nu, "GoodBeamDetPOTCounting");
    
    // Is the coil okay? Only cut if we have the cutOnDataQuality flag
    // This seems to be redundant with the GoodBeamDetPOT, but this is what
    // seems to have been implemented.
    Keep_Data_If(nu.coilIsOk || !nu.cutOnDataQuality, nu, "CoilIsOkay");

    // Filter out LI events
    Cut_If(nu.isLI || nu.litime != -1, "IsLI");    
  
    // Cut on good beam (IsGoodBeam)
    Keep_Data_If(CC0325Std::IsGoodBeam(nu) || !nu.cutOnBeamInfo, nu, "IsGoodBeam");
    
    // Per-detector cuts
    if (nu.detector == Detector::kFar)
    {
      // Data Quality
      Keep_Data_If(nu.isGoodDataQuality || !nu.cutOnDataQuality, nu, "DataQuality");
      
      // IsGoodTimeToNearestSpill
      // Use the improved JJE cut here
      Keep_Data_If(NuCuts::GoodTimeToNearestSpill(nu, -2, +12) || !nu.cutOnSpillTiming, nu, "SpillTime");
    
    }
    else if (nu.detector == Detector::kNear)
    {
      // Cut on the coil current
      Keep_Data_If(nu.coilCurrent < -1000 || !nu.cutOnDataQuality, nu, "CoilCurrent");
      
    } // End of per-detector cuts    
  }
  
  // The Main Selection ////////////////////////////////////////////////
  void CCAPresel::Selection(const NuEvent &)
  {
    // Since this is a preselection only class, we have no selection
  }

// A Base CC selection for the CC2010 analysis
  CCA::CCA(const NuPlots *plots, TString ana) : 
    NuCut("CCA", plots), rescuts(0), cutstring("")
  {
    SetFidVol("cc2008");
    cutstring=ana;
    SetCutFile("/minos/data/users/sjc/kNN_MicroDST/cuts_numu_D07d3.root");
    MSG("NuCutImps",Msg::kInfo)
      <<"Initializing resolution binning for "<<cutstring<<endl;
  }
  
  void CCA::SetCutFile(const char* filename)
  {
    TDirectory *tmpd = gDirectory; //get the right path
    TFile *cutfile = new TFile(filename,"READ");
    gDirectory = tmpd;
    if(cutfile->IsZombie()){
      MSG("NuCutImps",Msg::kWarning)<<"cutfile does not exist" << endl;
      return;
    }
    rescuts = (TTree *)cutfile->Get("cuts");
    MSG("NuCutImps",Msg::kWarning)<<"cutfile exists: "
                                  <<rescuts->GetEntries()<<" entries"<<endl;
  }
  
  // The CC selection uses a slightly altered definition of fiducial volume
  // origin: NuCuts::IsInFidVolTrk circa Sept '09
  Bool_t CCA::InFidVol(const NuEvent &nu) const
  {
    // Move trk vtx upstream by 3.92cm from scintillator to steel
    // for this analysis version
    NuEvent nuc = nu;
    nuc.zTrkVtx = nu.zTrkVtx - (0.0392*Munits::m);
    
    return NuCut::InFidVol(nuc);
  }
  
  // The Preselection //////////////////////////////////////////////////
  void CCA::Preselection(const NuEvent &nu)
  {
    // Do the base preselection (shared between CC and NC)
    Defer_Preselection(fCCAPresel, nu);
    
    // Now do the rest of the preselection
    
    // Good tracks
    Keep_If(nu.ntrk >= 1, "GoodNumberofTracks");
    
    // Is it in the fiducial volume?
    // ntrk cut to suppress warnings
    Keep_If(nu.ntrk > 0 && InFidVol(nu), "FidVol");

    // Per-detector cuts
    if (nu.detector == Detector::kFar)
    {
      // TrackFitPass
      Keep_If(nu.trkfitpass == 1, "trackfitpass");
      
      // Cut on track direction angle
      Keep_If(nu.dirCosNu > 0.6);
      
    }  
    else if (nu.detector == Detector::kNear)
    {
      // Need to implement ND Trackfitpass (uses reclamation)
      // See Zeynep talk DocDB-6382 - mentions 'Nikis Reclamation'
      Keep_If(NuCuts::IsGoodTrackFitPassReclamation(nu), "trackfitreclamation");
    }
  }
  
  // The Main Selection ////////////////////////////////////////////////
  void CCA::Selection(const NuEvent &nu)
  {
    Keep_If(nu.roID>0.3);
    
    if(cutstring=="CCA") return;
    
    if(cutstring=="CCA_NUBAR"){
      Keep_If(nu.charge>0);
    }
    
    else{
      if(nu.detector == Detector::kFar/* && nu.charge<0*/){
        Resolution(nu);
      }
    }
  }
  
  void CCA::Resolution(const NuEvent &nu)
  {
    Double_t low_en_edge=0.0;
    Double_t high_en_edge=0.0;
    Double_t low_res_edge=0.0;
    Double_t high_res_edge=0.0;
    Double_t quantiles[6];
    
    rescuts->SetBranchAddress("quantiles",quantiles);
    rescuts->SetBranchAddress("lowedge",&low_en_edge);
    rescuts->SetBranchAddress("highedge",&high_en_edge);
    
    rescuts->GetEntry(rescuts->GetEntries());
    

    for(Int_t i=0;i<rescuts->GetEntries();i++){
      rescuts->GetEntry(i);
      
      if(cutstring=="CCA_0"){
        low_res_edge=0.0;
        high_res_edge=quantiles[1];
      }
      if(cutstring=="CCA_1"){
        low_res_edge=quantiles[1];
        high_res_edge=quantiles[2];
      }
      if(cutstring=="CCA_2"){
        low_res_edge=quantiles[2];
        high_res_edge=quantiles[3];
      }
      if(cutstring=="CCA_3"){
        low_res_edge=quantiles[3];
        high_res_edge=quantiles[4];
      }
      if(cutstring=="CCA_4"){
        low_res_edge=quantiles[4];
        high_res_edge=9999999;
      }
      
      if(nu.energy>low_en_edge && nu.energy<high_en_edge){
        MAXMSG("NuCutImps",Msg::kInfo,20)
          <<"Entry:"<<i<<" Energy:"<<nu.energy<<" sigma/E:"<<(nu.resolution/nu.energy)<<endl
          <<"Keep if between "<<low_res_edge<<" and "<<high_res_edge<<endl;
        Keep_If((nu.resolution/nu.energy)>low_res_edge && (nu.resolution/nu.energy)<high_res_edge,"resolution");
        return;
      }
    }
  }

  // A Base NC selection for the CC2010 analysis

  CCA_NC::CCA_NC(const NuPlots* plots) : 
    NuCut("CCA_NC", plots)
  {
    // Match CC fiducial volume for simplicity
    // TODO CJB - will this try to use the track vertex? Is that bad?
    SetFidVol("cc2008");
  }
  
  // Fiducial Volume ///////////////////////////////////////////////////
  Bool_t CCA_NC::InFidVol(const NuEvent& nu) const
  {
    // If there's a track and it's longer than the shower then use its
    // vertex, otherwise use the event vertex. This is what is done in
    // the NC analysis.
    // TODO CJB - is this maximally consistent with the CC analysis?
    if(nu.ntrk > 0 && nu.trknplane > nu.nplaneShw){
      // Move trk vtx upstream by 3.92cm from scintillator to steel
      NuEvent nuc = nu;
      nuc.zTrkVtx = nu.zTrkVtx - (0.0392*Munits::m);
      return NuCut::InFidVol(nuc);
    }
    return NuCut::InFidVolEvt(nu);
  }

  // The Preselection //////////////////////////////////////////////////
  void CCA_NC::Preselection(const NuEvent& nu)
  {
    Defer_Preselection(fCCAPresel, nu);

    // Now do the rest of the preselection

    // No good track cut

    // Is it in the fiducial volume?
    Keep_If(InFidVol(nu), "FidVol");

    // No track fit pass or direction cosines or track reclamation

    // TODO CJB - check CC cuts don't overcut these
    if(nu.detector == Detector::kFar ) FDCleaning(nu);
    if(nu.detector == Detector::kNear) NDCleaning(nu);
  }
  
  // The Main Selection ////////////////////////////////////////////////
  void CCA_NC::Selection(const NuEvent& nu)
  {
    Keep_If(nu.roID < 0.3, "RoID");
  }
  
  // FD Cleaning ///////////////////////////////////////////////////////
  void CCA_NC::FDCleaning(const NuEvent& nu)
  {
    // See NCUtils/Cuts/NCAnalysisCutsNC for original version

    Cut_If(nu.nevt == 0, "NumEvents");

    const double phFrac = nu.evtphsigcor/max(nu.snarlPulseHeight, 1e-10);
      Cut_If(nu.nevt > 2 || (nu.nevt == 2 && phFrac < .75), "PHfrac");

    FDCosmics(nu); // Cut out cosmics

    Cut_If(nu.evtphsigcor < 5000 && nu.nstripEvt <= 4, "FibreNoise");

    // Timing cut. Data only
    const double dt_spill = nu.trigtime*Munits::s - nu.nearestSpillNanosec*Munits::ns;

    Cut_Data_If(dt_spill < -2*Munits::microsecond ||
                dt_spill > 12*Munits::microsecond,
                nu, "Timing");
  }

  template<class T> inline T SQR(T t){return t*t;}

  // FD Cleaning : cosmics ///////////////////////////////////////////////
  void CCA_NC::FDCosmics(const NuEvent& nu)
  {
    //Ultra-steep showers
    Cut_If(nu.planeEvtN == 0 || nu.nstripEvt >= SQR(nu.planeEvtN),
           "UltraSteep");

    //Steep showers
    const double tRMS = TMath::Sqrt(SQR(nu.transverseRMSU)+
                                    SQR(nu.transverseRMSV));
    Cut_If(nu.nshw > 0 && tRMS > .3+.1901*TMath::Log10(nu.planeEvtN), "Steep");

      //Track angle
    Cut_If(nu.ntrk > 0 && nu.trkvtxdcosz < 0.4, "TrkDCosZ");

      const bool upExitTrack = ((nu.zTrkEnd > 28.78 ||
                                 nu.endMetersToCloseEdge < .5) &&
                                nu.yTrkEnd > -1.657);

    Cut_If(nu.ntrk > 0 && upExitTrack && nu.dtdz < -.5 && nu.trkvtxdcosy < -.4,
             "UpExitTrack");
  }

  // ND Cleaning ///////////////////////////////////////////////////////
  void CCA_NC::NDCleaning(const NuEvent& nu)
  {
    // See NCUtils/Cuts/NCAnalysisCuts::IsMultiCutsClean for original version

    // 35ns timing cut
    Cut_If(TMath::Abs(nu.minTimeSeparation) < 35e-9, "TimeSep");

    // tiny events are mostly junk
    Cut_If(nu.nstripEvt < 5, "TotalStrips");

    // this cuts very steep showers, leaking in
    Cut_If(nu.nstripEvt > 1.15*SQR(nu.planeEvtN), "Steep");

    // this cuts leakage which leaves activity in partially instrumented
    // region
    Cut_If(nu.edgeActivityStrips > 2 &&
           nu.edgeActivityPH > 1000 &&
           nu.energy < 5 &&
           nu.nplaneShw > nu.trknplane,
           "Leakage");

    Cut_If(nu.oppEdgeStrips > 2 &&
           nu.energy < 5.0 &&
           nu.nplaneShw > nu.trknplane,
           "OppositeEdge");

    // make additional deltaZ cuts if (|minDeltaT|<175ns)
    Cut_If(TMath::Abs(nu.closeTimeDeltaZ) < 1.0 &&
           TMath::Abs(nu.minTimeSeparation) < 175e-9,
           "dz&dt");
  }

// NuBar front face rock event cross check by Matt Strait 2009-10
  ChairSound::ChairSound(const NuPlots * plots):
    NuCut("ChairSound", plots)
    // fNuCutsSelection(NuCuts::kNMB0325ChairSound, plots)
  {
    SetFidVol("cc2008");
  }

  // Recalculate the region of the detector!
  Int_t ChairSound::GetRegion(const NuEvent &nu)
  {
    // An enum for the rock detector region - this probably shouldn't go here,
    // but for now it will have to do (would rather use an enum than plain numbers)
    enum Rock_DetectorRegion {
      kNoTrack   = -1,
      kFiducial  = 0,
      kFront     = 1,
      kBack      = 2,
      kGap       = 3,
      kEdge      = 4,
      kCoilHole  = 5,
      kCorner    = 6,
      kUnknown   = 7
    };
  
    // Placeholder to avoid nesting if's in this function
    Bool_t inFid = false, nearFront = false, nearEdge = false;
    Bool_t nearGap = false, nearCoil = false, nearBack = false;
    Double_t radius2 = nu.xTrkVtx*nu.xTrkVtx + nu.yTrkVtx*nu.yTrkVtx;
    Double_t z = nu.zTrkVtx - FidVol::getTrkVtxZOffset();
    // Get the detector
    Detector::Detector_t det = (Detector::Detector_t)nu.detector;
  
    // If we are the near detector, we can't calculate the region - return unknown
    if (det == Detector::kNear) return kUnknown;

    // Now - examine each of the regions and determine if the vertex is in there
    // 0 = Fiducial
    inFid = infid((Detector::Detector_t)nu.detector, (SimFlag::SimFlag_t)nu.simFlag,
                  nu.xTrkVtx, nu.yTrkVtx, z);
    if (inFid) return kFiducial;
  
    // 1 = Front, non-edge
    nearFront = z < FidVol::getFarZMC(0);
    nearEdge = radius2 >= (FidVol::getFarRouter()*FidVol::getFarRouter());
    nearCoil = radius2 < FidVol::getFarRinner()*FidVol::getFarRinner();
  
    if(nearFront && !nearEdge && !nearCoil) return kFront;

    // 2 = Back, anywhere (no corners at back)
    nearBack = z > FidVol::getFarZMC(3);
    if(nearBack) return kBack;

    // 3 = Supermodule gap, non-edge
    nearGap = z > FidVol::getFarZMC(1) && z < FidVol::getFarZMC(2);
    if(nearGap && !nearEdge && !nearCoil) return kGap;

    // 4 = edge, non-end
    if(nearEdge && !nearBack && !nearFront && !nearGap) return kEdge;

    // 5 = coil hole, non-end
    if(nearCoil && !nearBack && !nearFront && !nearGap) return kCoilHole;

    // 6 = corner cases
    // No other cases have been matched, it must (?) be a corner case?
    return kCorner;
  }
  
  void ChairSound::Preselection(const NuEvent &nu)
  {
    // can't use CC0325Std::Preselection because it has a fiducial cut and
    // an overly tight angle cut.  This is copied from there except for those.

    Keep_If(nu.ntrk >= 1, "GoodNumberofTracks");    // Good tracks
    Cut_If(nu.isLI || nu.litime != -1, "LICut");    // Filter out LI events
    Keep_Data_If(!nu.cutOnDataQuality || nu.coilIsOk, nu, "IsCoilOkay");
    Keep_Data_If(nu.goodBeamSntp && nu.coilIsOk, nu,"GoodBeamDetPOT");
    Keep_Data_If(!nu.cutOnBeamInfo || CC0325Std::IsGoodBeam(nu),nu,"GoodBeam");

    // Keep only front face events that are neither near the edge nor coil hole
    // Keep_If(nu.regionTrkVtx == 1, "front_face");
    Keep_If(GetRegion(nu) == 1, "front_face");

    // Per-detector cuts
    if (nu.detector == Detector::kFar){
      Keep_If(nu.trkfitpass == 1, "trackfitpass");
      Keep_Data_If(nu.isGoodDataQuality||!nu.cutOnDataQuality,nu,"DataQuality");

      // [-2, 12] microseconds
      Keep_Data_If(NuCuts::GoodTimeToNearestSpill(nu, -2, +12) 
                   || !nu.cutOnSpillTiming, nu, "SpillTime");
      // Cosmic cut?
      Keep_If(TMath::Abs(nu.trkvtxdcosy) < 0.4, "CosY");
    }
    else if (nu.detector == Detector::kNear){
      MAXMSG("CC0325Std",Msg::kWarning, 1) << "Front face numubar selection "
         " called on near detector. I doubt that this works properly.\n";
      Keep_If(NuCuts::IsGoodTrackFitPassReclamation(nu), "trackfitreclamation");
      Keep_Data_If(nu.coilCurrent<-1000||!nu.cutOnDataQuality,nu,"CoilCurrent");
    }
  }

  // The Main Selection ////////////////////////////////////////////////
  void ChairSound::Selection(const NuEvent &nu)
  {
    // Other selectors have a Defer_Selection to CC0325Std here, but I 
    // don't think it makes sense here since we don't share a preselection.
    // However, I honestly have no clear idea what it's even supposed to do, 
    // so I could be wrong.

    Keep_If(nu.roID >= 0.25, "RoID"); // Note: >= is not the same as >
    Keep_If(1/nu.sigqp_qp > 3.5, "SigmaQP_QP");
    Keep_If(TMath::Abs(nu.relativeAngle - TMath::Pi()) > 2.08, "RelativeAngle");
  }

// NuBar front face rock event cross check by Matt Strait 2009-10
  MSRock_Nov09::MSRock_Nov09(const NuPlots * plots):
    NuCut("MSRock_Nov09", plots)
  {
    SetFidVol("cc2008");
  }

  void MSRock_Nov09::Preselection(const NuEvent &nu)
  {
    // can't use CC0325Std::Preselection because it has a fiducial cut and
    // an overly tight angle cut.  This is copied from there except for those.

    Keep_If(nu.ntrk >= 1, "GoodNumberofTracks");    // Good tracks
    Cut_If(nu.isLI || nu.litime != -1, "LICut");    // Filter out LI events
    Keep_Data_If(!nu.cutOnDataQuality || nu.coilIsOk, nu, "IsCoilOkay");
    Keep_Data_If(nu.goodBeamSntp && nu.coilIsOk, nu,"GoodBeamDetPOT");
    Keep_Data_If(!nu.cutOnBeamInfo || CC0325Std::IsGoodBeam(nu),nu,"GoodBeam");

    // Per-detector cuts
    if (nu.detector == Detector::kFar){
      // Matt's rock selection does not use Trackfitpass at the moment
      // Keep_If(nu.trkfitpass == 1, "trackfitpass");
      Keep_Data_If(nu.isGoodDataQuality||!nu.cutOnDataQuality,nu,"DataQuality");

      // [-2, 12] microseconds
      Keep_Data_If(NuCuts::GoodTimeToNearestSpill(nu, -2, +12) 
                   || !nu.cutOnSpillTiming, nu, "SpillTime");
      // Cosmic cut?
      Keep_If(TMath::Abs(nu.trkvtxdcosy) < 0.4, "CosY");
    }
    else if (nu.detector == Detector::kNear){
      MAXMSG("CC0325Std",Msg::kWarning, 1)
        << "Rock Muon Selection being called on ND. You probably don't want to be doing this."
        << endl;
      Keep_If(NuCuts::IsGoodTrackFitPassReclamation(nu), "trackfitreclamation");
      Keep_Data_If(nu.coilCurrent<-1000||!nu.cutOnDataQuality,nu,"CoilCurrent");
    }
    
    // Make a fiducial cut
    // Keep_If(nu.regionTrkVtx == 1, "front_face");
    // Keep_If(GetRegion(nu) == 1, "front_face");
    Keep_If(ChairSound::GetRegion(nu) > 0 && ChairSound::GetRegion(nu) < 7, "Antifiducial");
  }

  // The Main Selection ////////////////////////////////////////////////
  void MSRock_Nov09::Selection(const NuEvent &nu)
  {
    // Other selectors have a Defer_Selection to CC0325Std here, but I 
    // don't think it makes sense here since we don't share a preselection.
    // However, I honestly have no clear idea what it's even supposed to do, 
    // so I could be wrong.
    // Keep_If(nu.roID >= 0.25, "RoID"); // Note: >= is not the same as >
    // Keep_If(1/nu.sigqp_qp > 3.5, "SigmaQP_QP");
    // Keep_If(TMath::Abs(nu.relativeAngle - TMath::Pi()) > 2.08, "RelativeAngle");

    Keep_If(nu.roID >= 0.1, "RoID");
  }




} // Closes the namespace

/*


*/

---