NueStandard.cxx

Go to the documentation of this file.

// $Id: NueStandard.cxx,v 1.59 2009/10/27 16:19:05 scavan Exp $
//
// The NueConvention File will hopefully organize
//  frequently used nue standards such as
//    event class (was ClassType)
//    fiducial Volume
//    file indexing
//
// Author: Josh Boehm
// Created: Sept. 14, 2006
#include <vector>
#include "NueAna/NueAnaTools/NueConvention.h"
#include "NueStandard.h"
#include "DcsUser/CoilTools.h"
#include "MCReweight/SKZPWeightCalculator.h"
#include "TMultiLayerPerceptron.h"

#include "TFile.h"
#include "TMath.h"
#include <cmath>

#include <fstream>
#include "MessageService/MsgService.h"

CVSID("");

bool NueStandard::IsInFid(NueRecord *nr)
{
   float x = nr->srevent.vtxX;
   float y = nr->srevent.vtxY;
   float z = nr->srevent.vtxZ;
                                                                               
   Detector::Detector_t fDet;
   fDet = nr->GetHeader().GetVldContext().GetDetector();
   if (fDet != Detector::kFar && fDet != Detector::kNear) return false;

   SimFlag::SimFlag_t sim = nr->GetHeader().GetVldContext().GetSimFlag();
   bool isMC = (sim == SimFlag::kMC);
   int retval = 0;

   if (fDet == Detector::kFar)
     retval = NueConvention::IsInsideFarFiducial_Nue_Standard(x,y,z, isMC);
 
   if (fDet == Detector::kNear)
     retval = NueConvention::IsInsideNearFiducial_Nue_Standard(x,y,z, isMC);
      
   return (retval == 1);
}

bool NueStandard::PassesPOTStandards(NueRecord *nr)
{
   SimFlag::SimFlag_t sim = nr->GetHeader().GetVldContext().GetSimFlag();
   bool isMC = (sim == SimFlag::kMC);
   VldContext vld = nr->GetHeader().GetVldContext();

   if (isMC) return true;
   
   //common to both
   if (nr->bmon.goodBeamMon != 1) return false;

   if (nr->GetHeader().GetVldContext().GetDetector() == Detector::kNear )
   {
     bool goodCoil =  nr->eventq.coilQuality && (nr->eventq.coilDirection ==1);
                      //CoilTools::IsOK(vld) && !CoilTools::IsReverse(vld);
     int dpDQ = nr->eventq.passNearDetQuality;
     return goodCoil && (dpDQ == 1);
   }

   if (nr->GetHeader().GetVldContext().GetDetector() == Detector::kFar)
   {

    if (nr->eventq.spillType != 1) return false;
    int dpDQ = nr->eventq.passFarDetQuality;
    if (nr->GetHeader().GetEventNo() < 0) return (dpDQ == 1);

     return (dpDQ == 1) && (nr->eventq.passLI == 1);
   }

   return false;
}                                                                                                          
bool NueStandard::PassesDataQuality(NueRecord *nr)
{                                                 
   SimFlag::SimFlag_t sim = nr->GetHeader().GetVldContext().GetSimFlag();
   bool isMC = (sim == SimFlag::kMC);
   VldContext vld = nr->GetHeader().GetVldContext();

   if (isMC) return true;
   
   if(!PassesPOTStandards(nr)) return false;

   if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kNear) return true;

   if (nr->GetHeader().GetVldContext().GetDetector() == Detector::kFar){

     int rcB = nr->eventq.rcBoundary;
     float tmin = nr->srevent.eventTimeMin;
     double spillT = nr->bmon.dt_stnd;   
     return PassesFarDataQuality(-1, rcB, 1, tmin, spillT);       
   }

   return false;
}

bool NueStandard::PassesNearDataQuality(int gbm, float cc, int st)
{
  //Assumes it is taking:
  // nr->bmon.goodBeamMon, nr->srevent.coilCurrent, nr->srevent.spillType
  if (gbm == 1 && cc < -1000.0 && st != 3)
     return true;
                                                                                
  return false;
}

bool NueStandard::PassesFarDataTiming(NueRecord *nr) 
{
 if(
       Detector::kFar == nr->GetHeader().GetVldContext().GetDetector()
    && SimFlag::kData == nr->GetHeader().GetVldContext().GetSimFlag()
  )
 {
  float tmin = nr->srevent.eventTimeMin;
  double spillT = nr->bmon.dt_stnd;
  return( (tmin - spillT) > -2e-6 && (tmin - spillT) < 12e-6 );
 }
  return true;
}

bool NueStandard::PassesFarDataQuality(float li, int rc, int dpfddq, 
                                          float tmin, double spillt)
{
  if (li == -1 && rc == 0 && dpfddq == 1
    && (tmin - spillt) > -20e-6 && (tmin - spillt) < 30e-6)
     return true;
                                                                                
  return false;
}

bool NueStandard::IsGoodRun(NueRecord *nr)
{
  if(nr->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)
  {
    return true;
  }
  
  if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kFar)
  {
    return NueStandard::IsGoodFarRun(nr);
  }
  
  return false;
}

bool NueStandard::IsGoodNearRun(NueRecord *nr)
{
  if(nr->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)
  {
    return true;
  }
  
  if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kFar)
  {
    return true;
  }
  
  //for near detector data
  if(nr->GetHeader().GetRun()==8165 || nr->GetHeader().GetRun()==11318)
  {
    return false;
  }
  
  if(nr->GetHeader().GetRun()==7942 && nr->GetHeader().GetSubRun()==0)
  {
    return false;
  }
  
  if(nr->GetHeader().GetRun()==7982 && nr->GetHeader().GetSubRun()==0)
  {
    return false;
  }
  
  if(nr->GetHeader().GetRun()==8214 && nr->GetHeader().GetSubRun()==0)
  {
    return false;
  }
  
  if(nr->GetHeader().GetRun()==9809 && nr->GetHeader().GetSubRun()==0)
  {
    return false;
  }
  
  return true;
}

bool NueStandard::IsGoodFarRun(NueRecord *nr)
{
  if(nr->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)
  {
    return true;
  }
  
  if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kNear)
  {
    return true;
  }
  
  //for far detector data
  static vector<FilePosition> badFDrunlist;
  static unsigned int listpos = 0;
  int Run, SubRun;
  
  // If this is the first call fill the file list
  if(badFDrunlist.size() == 0)
  {
    std::ifstream ins;
    char *srt_dir = getenv("SRT_PRIVATE_CONTEXT");
    ins.open(Form("%s/NueAna/badrunsFD_list.txt",srt_dir));
    if(!ins.is_open())
    {
      MSG("NueStandard", Msg::kError)<<"IsGoodFarRun(): Unable to open badrunsFD_list.txt"<<endl;
    }

    while(!ins.eof())
    {
      ins>>Run>>SubRun;
      FilePosition temp;
      if(!ins.eof()){
          temp.Run = Run;
          temp.SubRun = SubRun;
          temp.Snarl = -1;
          temp.Event = -1;
          badFDrunlist.push_back(temp);
      }
    }
    
    if(badFDrunlist.size() == 0)  return true;
   }
   
   FilePosition current;
   
   current.Snarl = -1;
   current.Event = -1;
   current.Run = nr->GetHeader().GetRun();
   current.SubRun = nr->GetHeader().GetSubRun();
   
   if(current < badFDrunlist[0]) return true;

   if(current > badFDrunlist[badFDrunlist.size() - 1]) return true;

   if(current > badFDrunlist[listpos] )
   {
      while(current > badFDrunlist[listpos] && listpos < badFDrunlist.size())
          listpos++;
   }
   else if(current < badFDrunlist[listpos] )
   {
      while(current < badFDrunlist[listpos] && listpos > 0)
          listpos--;
   }

   if(current == badFDrunlist[listpos])
   {
       return false;
   }

   
   return true;
}

//  Now to the PreSelection Cuts

bool NueStandard::PassesTrackPlaneCut(NueRecord *nr)
{
   int tp = TMath::Abs(nr->srtrack.endPlane - nr->srtrack.begPlane);

   return PassesTrackPlaneCut(tp);
}

bool NueStandard::PassesTrackPlaneCut(int trkplane)
{
   return (trkplane < 25);
}

bool NueStandard::PassesMinPlaneCut(NueRecord *nr)
{
   int planes = nr->shwfit.contPlaneCount050;
   return PassesMinPlaneCut(planes);
}

bool NueStandard::PassesMinPlaneCut(int planes)
{
   return (planes > 4);
}

bool NueStandard::PassesShowerCut(NueRecord *nr)
{
  int nshw = nr->srevent.showers;
  return PassesShowerCut(nshw);
}

bool NueStandard::PassesShowerCut(int nshw)
{
  return (nshw > 0);
}

bool NueStandard::PassesCosmicCutFunction(NueRecord *nr)
{
   bool result=true;    
   //cosmicy shower cut
   result = result && !(nr->shwfit.UVSlope > 10);

   //cosmicy track cut
   if (nr->srevent.tracks<1) return result;  //require a track to apply this cut

   float ex =  nr->srtrack.endX;
   float ey =  nr->srtrack.endY;
   float ez =  nr->srtrack.endZ;

   float vx =  nr->srtrack.vtxX;
   float vy =  nr->srtrack.vtxY;
   float vz =  nr->srtrack.vtxZ;

   float cosy = acos(fabs(vy-ey)/sqrt((vx-ex)*(vx-ex)+(vy-ey)*(vy-ey)+(vz-ez)*(vz-ez))); 
   float dy=fabs(vy-ey);

   result = result && !(dy>2 && cosy < 0.6);

   return result;
}

bool NueStandard::PassesCosmicCut(NueRecord *nr)
{
  if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kNear)
       return true;

   bool ret=true;
   if(nr->eventq.passCosmicCut > -10)
     return (nr->eventq.passCosmicCut == 1);

   //or we are still in Ent
 
   if(nr->dtree.bt_var1==0)ret = false;
   else if(nr->dtree.bt_var1==1)ret = true;
   else{ NueStandard::FillCosmicCut(nr); ret = nr->dtree.bt_var1; }

   Detector::Detector_t fDet;
   fDet = nr->GetHeader().GetVldContext().GetDetector();
   if (fDet != Detector::kFar) ret = true;

   return ret;
}

void NueStandard::FillCosmicCut(NueRecord *nr)
{
   int bv = 0;
   if (NueStandard::PassesCosmicCutFunction(nr)) bv=1;
   nr->dtree.bt_var1=bv;
}

bool NueStandard::PassesTrackLikePlaneCut(NueRecord *nr)
{
   int tlp = nr->srtrack.trklikePlanes;
   return PassesTrackLikePlaneCut(tlp);
}

bool NueStandard::PassesTrackLikePlaneCut(int tlp)
{       
   return (tlp < 16);
}

bool NueStandard::PassesLowEnergyCut(NueRecord *nr)
{
   float energy = nr->srevent.phNueGeV;
   return PassesLowEnergyCut(energy);
}
 
bool NueStandard::PassesLowEnergyCut(float energy)
{
   return (energy > 1.0);
}

bool NueStandard::PassesHighEnergyCut(NueRecord *nr)
{
   float energy = nr->srevent.phNueGeV;
   return PassesHighEnergyCut(energy);
}
          
bool NueStandard::PassesHighEnergyCut(float energy)
{
   return (energy < 8.0);
}

// Now to combinations of Preselection Cuts                                                                                                          
bool NueStandard::PassesPreSelection(NueRecord *nr)
{
   bool temp = PassesNonHEPreSelection(nr);
   temp = temp && PassesHighEnergyCut(nr);
  
   return temp;
}

bool NueStandard::PassesPreSelection(int trkplane, int tlp, float energy)
{
   bool temp = PassesNonHEPreSelection(trkplane, tlp, energy);
   temp = temp && PassesHighEnergyCut(energy);
   return temp;
}

bool NueStandard::PassesNonHEPreSelection(NueRecord *nr)
{
   bool temp = PassesPreSelectionTrackCuts(nr);
   temp = temp && PassesLowEnergyCut(nr);
   temp = temp && PassesPreSelectionBasicCuts(nr);
   
   return temp;
}                                                                                                       
bool NueStandard::PassesNonHEPreSelection(int trkplane, int tlp, float energy)
{
   bool temp = PassesPreSelectionTrackCuts(trkplane, tlp);
   temp = temp && PassesLowEnergyCut(energy);
   return temp;
}

bool NueStandard::PassesPreSelectionTrackCuts(NueRecord *nr)
{
   bool temp = PassesTrackPlaneCut(nr);
   temp = temp && PassesTrackLikePlaneCut(nr);
   return temp;  
}

bool NueStandard::PassesPreSelectionTrackCuts(int trkplane, int tlp)
{
   bool temp = PassesTrackPlaneCut(trkplane);
   temp = temp && PassesTrackLikePlaneCut(tlp);
   return temp;
}

bool NueStandard::PassesPreSelectionBasicCuts(NueRecord *nr)
{
   bool one = PassesMinPlaneCut(nr);
   one = one && PassesCosmicCut(nr); 
   one = one && IsLargestEvent(nr);
   one = one && PassesShowerCut(nr);

   return one;
}

bool NueStandard::PassesMRCCFiducial(NueRecord *nr)
{
   float x = nr->mri.vtxx;
   float y = nr->mri.vtxy;
   float z = nr->mri.vtxz;
                                                                                
   Detector::Detector_t fDet;
   fDet = nr->GetHeader().GetVldContext().GetDetector();
   
   SimFlag::SimFlag_t sim = nr->GetHeader().GetVldContext().GetSimFlag();
   bool isMC = (sim == SimFlag::kMC);
   int retval = 0;
                                                                                
   if (fDet == Detector::kFar)
   {
     if(sqrt((x*x) + (y*y))>0.3 && sqrt((x*x) + (y*y))<sqrt(15.5) && ((z>0.3 && z<14.4)||(z>16.1 && z<28.2)))
     {
       retval=1;
     }
   }
       
   if (fDet == Detector::kNear)
     retval = NueConvention::IsInsideNearFiducial_MRE_Standard(x,y,z, isMC);

   return (retval == 1);
}

bool NueStandard::PassesMRCCPreSelection(NueRecord *nr)
{
    double pid = nr->mri.orig_roCCPID;
    int fitp = nr->mri.fitp;
    float bestC = nr->mri.best_complete;

    return NueStandard::PassesMRCCPreSelection(bestC, fitp, pid);
}

bool NueStandard::PassesMRCCPreSelection(float bestComp, int fitp, double pid)
{
   //currently assumes RO-PID
   bool good = true;
   good = good && (bestComp > -10);
   good = good && (fitp == 1);
   good = good && (pid > 0.3);
                                                                                
   return good;
}

bool NueStandard::PassesMREFiducial(NueRecord *nr)
{
    return NueStandard::PassesMRCCFiducial(nr);
}

bool NueStandard::PassesMREPreSelection(NueRecord *nr)
{
    return NueStandard::PassesMRCCPreSelection(nr);
}
                                                                                                              
bool NueStandard::PassesMREPreSelection(float bestC, int fitp, double pid)
{
    return NueStandard::PassesMRCCPreSelection(bestC, fitp, pid);
}
                                                                                             
bool NueStandard::PassesSysPreSelectionNoHE(NueRecord *nr)
{
   int tp = TMath::Abs(nr->srtrack.endPlane - nr->srtrack.begPlane);
   int tlp = nr->srtrack.trklikePlanes;
   float energy = nr->srevent.phNueGeV;

   return PassesSysPreSelectionNoHE(tp, tlp, energy);
}
                                                                                             
bool NueStandard::PassesSysPreSelectionNoHE(int trkplane, int tlp, float energy)
{
   bool good = true;
   good = good && (trkplane < 28);
   good = good && (tlp < 18);
   good = good && (energy > 0.5);
   return good;
}


bool NueStandard::PassesSysPreSelection(NueRecord *nr)
{
   int tp = TMath::Abs(nr->srtrack.endPlane - nr->srtrack.begPlane);
   int tlp = nr->srtrack.trklikePlanes;
   float energy = nr->srevent.phNueGeV;
                                                                                
   return PassesSysPreSelection(tp, tlp, energy);
}
                                                                                
bool NueStandard::PassesSysPreSelection(int trkplane, int tlp, float energy)
{
   bool good = true;
   good = good && (trkplane < 28);
   good = good && (tlp < 18);
   good = good && (energy > 0.5 && energy < 10);
   return good;
}

double NueStandard::GetPIDValue(NueRecord *nr, Selection::Selection_t sel)
{
  double retval = -9999;

  switch(sel) {
  case Selection::kCuts:  retval = nr->treepid.fCutPID;           break;
  case Selection::kANN6:  retval = nr->ann.pid_6inp;              break;
  case Selection::kANN30:  retval = nr->ann.pid_30inp;            break;
  case Selection::kANN2PE:  retval = nr->ann.pid_11inp;            break;
  case Selection::kANN2PE_DAIKON04:  retval = nr->ann.pid_11inp_daikon04;         break;
  case Selection::kANN14_DAIKON04:   retval = nr->ann.pid;         break;
//  case Selection::kANN14_DAIKON04:   retval = nr->ann.pid_14inp_daikon04;         break;
  case Selection::kSSPID: retval = nr->subshowervars.pid;  break;
  case Selection::kMCNN:
     if (nr->mcnnv.bestmatches < 50)  retval = -0.5;
     else        retval = nr->mcnnv.mcnn_var1;                        break;   
  
  
  case Selection::kParticlePID:
        if(PassesParticlePIDPreselectionCut(nr))retval = nr->precord.event.pidF;
        break;  
  
  default:                                                        break;
  }

  return retval;
}

bool NueStandard::PassesPIDSelection(NueRecord *nr, Selection::Selection_t sel)
{
  bool pass = false;
  double val = GetPIDValue(nr, sel);
                                                                                                  
  switch(sel) {
  case Selection::kCuts:  pass = (int(val) == 1);       break;
  case Selection::kANN6:   pass = val > 0.5;            break;
  case Selection::kANN30:  pass = val > 0.5;            break;
  case Selection::kANN2PE:  pass = val > 0.7;            break;
  case Selection::kANN2PE_DAIKON04:  pass = val > 0.7;   break;
  case Selection::kANN14_DAIKON04:   pass = val > 0.75;  break;
  case Selection::kSSPID: pass = val > 0.67;             break;
  case Selection::kMCNN:  pass = val > 0.80;             break;
  case Selection::kParticlePID:  pass = PassesParticlePIDCut(nr); break;
  default:                                              break;
  }
                                                                                                  
  return pass;
}

bool NueStandard::PassesParticlePIDCut(NueRecord *nr)
{
        bool passPID=true;
        
        passPID = passPID && nr->precord.event.pidF>0.7;

        if(!passPID)return false;


        return passPID && PassesParticlePIDPreselectionCut(nr);
        

}


bool NueStandard::PassesParticlePIDPreselectionCut(NueRecord *nr)
{
                
        //fiducial cut
        bool passFid=true;
        passFid = passFid && (nr->precord.event.inFiducial==1); 

        if(!passFid)return false;


        double visenergy = nr->precord.event.visenergy;
                
        //preselection cut
        int det = nr->GetHeader().GetVldContext().GetDetector();
        if(det == Detector::kFar)
        {
                float offset =   0.489987;
                float slope  =  0.0387301;
        visenergy = visenergy*slope+offset;
        }else if(det == Detector::kNear)
        {
                float offset = 0.4803261;
                float slope = 0.03799819;
                visenergy = visenergy*slope+offset;
        }else{
                printf("please don't run in a detector that is not near or far (NueStandard::PassesParticlePIDCut)\n");
                exit(1);
        }
                
        bool passPre=true;
                
        double event_length=nr->precord.event.max_z-nr->precord.event.min_z;    
        passPre = passPre && event_length>0.1 && event_length<1.2;
        passPre = passPre && nr->precord.particles.longest_s_particle_s>0.1 && 
                nr->precord.particles.longest_s_particle_s<1.2;
        passPre = passPre && nr->precord.particles.ntot>0;
        passPre = passPre && visenergy>0.5 && visenergy<8;      

        return passPre && passFid;

}


bool NueStandard::IsLargestEvent(NueRecord *nr)
{
 if(nr->GetHeader().GetVldContext().GetDetector() == Detector::kFar){
    if(nr->srevent.largestEvent == 1) return true;
    return false;
 }
 return true;
}

bool NueStandard::PassesSelection(NueRecord *nr, Selection::Selection_t sel)
{
  bool pass = false;

  bool dq = NueStandard::PassesDataQuality(nr);
  bool infid = NueStandard::IsInFid(nr) && dq;
  
  bool isMR = nr->GetHeader().FoundMR();

  if (isMR){
    bool mrfid = PassesMRCCFiducial(nr);
    bool mrps = PassesMRCCPreSelection(nr);

    infid = infid && mrfid && mrps;
  } 
  
  bool basic = NueStandard::PassesPreSelectionBasicCuts(nr) && infid;
  bool presel = infid && NueStandard::PassesPreSelection(nr);
  bool pid = NueStandard::PassesPIDSelection(nr,sel); 

  switch(sel) {
  case Selection::kNone:     pass = true;                 break;
  case Selection::kDataQual: pass = dq;                   break;
  case Selection::kFid:      pass = infid;                   break;
  case Selection::kBasic:    pass = basic;               break;
  case Selection::kPre:   pass = presel;                  break;
  case Selection::kCuts:  pass = presel && pid;           break;
  case Selection::kANN6:  pass = presel && pid;          break;
  case Selection::kANN30: pass = presel && pid;          break;
  case Selection::kANN2PE: pass = presel && pid;          break;
  case Selection::kANN2PE_DAIKON04: pass = presel && pid; break;
  case Selection::kANN14_DAIKON04: pass = presel && pid;  break;
  case Selection::kSSPID: pass = presel && pid;           break;
  case Selection::kMCNN:  pass = presel && pid;           break;
  case Selection::kMDA:   pass = presel && false;         break;
  case Selection::kBDT:   pass = presel && false;         break;
  case Selection::kKNue:  pass = presel && false;         break;
  case Selection::kCC:    pass = infid && 
                     NueStandard::PassesCCSelection(nr);  break;
  default:                                                break;
  }

  return pass;
}

bool NueStandard::PassesSelection(NueRecord *nr, Selection::Selection_t sel, Selection::Selection_t sel2)
{
  bool pass = false;

  bool dq = NueStandard::PassesDataQuality(nr);
  bool infid = NueStandard::IsInFid(nr) && dq;
  bool basic = NueStandard::PassesPreSelectionBasicCuts(nr);
  bool presel = NueStandard::PassesPreSelection(nr);
  bool pid = NueStandard::PassesPIDSelection(nr,sel);

  switch(sel2) {
  case Selection::kDataQual:  dq = true;                  break;
  case Selection::kFid:       infid = true;               break;
  case Selection::kBasic:     basic = true;               break;
  case Selection::kPre:       presel = true;              break;
  default:                                                break;
  }
  
  bool isMR = nr->GetHeader().FoundMR();

  if(isMR){
    bool mrfid = true; //PassesMRCCFiducial(nr);
    bool mrps = PassesMRCCPreSelection(nr);

    infid = infid && mrfid && mrps;
  }
  
  basic = basic && infid;
  presel = dq && infid && presel;

  switch(sel) {
  case Selection::kNone:  pass = true;                    break;
  case Selection::kDataQual:  pass = dq;                  break;
  case Selection::kFid:   pass = infid;                   break;
  case Selection::kBasic: pass = basic;                   break;
  case Selection::kPre:   pass = presel;                  break;
  case Selection::kCuts:  pass = presel && pid;           break;
  case Selection::kANN6:   pass = presel && pid;          break;
  case Selection::kANN30:  pass = presel && pid;          break;
  case Selection::kANN2PE: pass = presel && pid;          break;
  case Selection::kANN2PE_DAIKON04: pass = presel && pid; break;
  case Selection::kANN14_DAIKON04: pass = presel && pid;  break;
  case Selection::kSSPID: pass = presel && pid;           break;
  case Selection::kMCNN:  pass = presel && pid;           break;
  case Selection::kMDA:   pass = presel && false;         break;
  case Selection::kBDT:   pass = presel && false;         break;
  case Selection::kKNue:  pass = presel && false;         break;
  case Selection::kCC:    pass = infid &&
                     NueStandard::PassesCCSelection(nr);  break;
  default:                                                break;
  }

  return pass;
}


bool NueStandard::PassesCCSelection(NueRecord *nr)
{
   int ntrack  = nr->srevent.tracks;
   int pass = nr->srtrack.passedFit;
   int endPlaneU = nr->srtrack.endPlaneU;
   int endPlaneV = nr->srtrack.endPlaneV;
   int deltaUVVtx = nr->srtrack.deltaUVVtx;   

   if (ntrack < 1){ pass = 0; endPlaneU = endPlaneV = deltaUVVtx = 300; }

   bool trackPass =  (pass == 1) || 
                (deltaUVVtx <= 5 && TMath::Abs(endPlaneU - endPlaneV) <= 40
                   && endPlaneU < 270 && endPlaneV < 270);


   return  (ntrack > 0) && trackPass && nr->anainfo.roCCPID > 0.3;
}


double NueStandard::GetRPWBeamWeight(NueRecord *nr, bool ismrcc)
{


  std::vector<double> pots;
  
  Detector::Detector_t det = nr->GetHeader().GetVldContext().GetDetector();
  BeamType::BeamType_t beam = nr->GetHeader().GetBeamType();

  if(beam != BeamType::kL010z185i
     &&beam != BeamType::kL010z185i_i124 
     &&beam != BeamType::kL010z185i_i191   
     &&beam != BeamType::kL010z185i_i213
     &&beam != BeamType::kL010z185i_i224
     &&beam != BeamType::kL010z185i_i232
     &&beam != BeamType::kL010z185i_i243
     &&beam != BeamType::kL010z185i_i257
     &&beam != BeamType::kL010z185i_i282
     &&beam != BeamType::kL010z185i_i303
     &&beam != BeamType::kL010z185i_i324
   
     &&beam != BeamType::kL010z000i
     &&beam != BeamType::kL010z000i_i209
     &&beam != BeamType::kL010z000i_i225
     &&beam != BeamType::kL010z000i_i232
     &&beam != BeamType::kL010z000i_i259
     &&beam != BeamType::kL010z000i_i300
     &&beam != BeamType::kL010z000i_i317
     &&beam != BeamType::kL010z000i_i326
     &&beam != BeamType::kL010z000i_i380){

    MAXMSG("NueStandard",Msg::kWarning,5)<<"Can't use NueStandard::GetRPWBeamWeight for non L010185 or L010000 beams, stop it!"<<endl;
    return 0;
  }

  if (det==Detector::kNear){
    if (ismrcc){
      for(int i=0;i<NRUNPERIODS;i++){
        pots.push_back(pot_ndmrcc[i]);
      }
    }
    else if (beam == BeamType::kL010z185i
     ||beam == BeamType::kL010z185i_i124
     ||beam == BeamType::kL010z185i_i191
     ||beam == BeamType::kL010z185i_i213
     ||beam == BeamType::kL010z185i_i224
     ||beam == BeamType::kL010z185i_i232
     ||beam == BeamType::kL010z185i_i243
     ||beam == BeamType::kL010z185i_i257
     ||beam == BeamType::kL010z185i_i282
     ||beam == BeamType::kL010z185i_i303
     ||beam == BeamType::kL010z185i_i324){
      for(int i=0;i<NRUNPERIODS;i++){
        pots.push_back(pot_nd[i]);
      }
    }
    else if (beam==BeamType::kL010z000i
     ||beam == BeamType::kL010z000i_i209
     ||beam == BeamType::kL010z000i_i225
     ||beam == BeamType::kL010z000i_i232
     ||beam == BeamType::kL010z000i_i259
     ||beam == BeamType::kL010z000i_i300
     ||beam == BeamType::kL010z000i_i317  
     ||beam == BeamType::kL010z000i_i326
     ||beam == BeamType::kL010z000i_i380){
      for(int i=0;i<NRUNPERIODS;i++){
        pots.push_back(pot_ndhornoff[i]);
      }
    }
  }
  else{
    if (ismrcc){
      for(int i=0;i<NRUNPERIODS;i++){
        pots.push_back(pot_fdmrcc[i]);
      }
    }
    //FD has no intensity samples
    else if (beam==BeamType::kL010z185i){
      for(int i=0;i<NRUNPERIODS;i++){
        pots.push_back(pot_fd[i]);
      }
    }
    //FD has no intensity samples
    else if (beam==BeamType::kL010z000i){
      MAXMSG("NueStandard",Msg::kWarning,5)<<"Can't use NueStandard::GetRPWBeamWeight for FD L010000 beams, stop it!"<<endl;
      return 0;
    }
  }

  if(nr->mctrue.nuFlavor < -1000){
    MAXMSG("NueStandard",Msg::kWarning,5)<<"No truth info for this event returning 0 for beamweight"<<endl;
    return 0;
  }

  double weight = GetRPWBeamWeight(nr->fluxweights.RPtotbeamweight,pots);
  return weight;
}

double NueStandard::GetRPWBeamWeight(std::vector<double> weights, std::vector<double> pots)
{

  double w=0.;
  double p=0.;
  for(unsigned int i=0;i<weights.size();i++){
    w+=weights[i]*pots[i];
    p+=pots[i];
  }

  if (p==0){
    cout<<"Total pots is set to zero, can not compute a exposure weighted average beam weight.  Returning 1"<<endl;
    return 1;
  }

  return w/p;
}

bool NueStandard::PassesNCCleaningCuts(NueRecord* nr)
{                                                                                
  // 40ns timing cut
  if (TMath::Abs(nr->eventq.minTimeSeparation)<40e-9)
    return false;
                                                                                
  // tiny events are mostly junk
  if (nr->srevent.totalStrips<5)    return false;
                                                                                
  // this cuts very steep showers, leaking in
  if ( (nr->srevent.totalStrips/(nr->srevent.planes*nr->srevent.planes))>1.0)
    return false;
                                                                                
  // this cuts leakage which leaves activity in partially instrumented
  // region
  if ( (nr->eventq.edgeActivityStrips>3)
       && (nr->eventq.edgeActivityPH>1000)
       && (nr->srevent.energyGeV<5.0)
       && (nr->srshower.planes>nr->srtrack.planes) ) // don't cut out CC
    return false;
  // also cut out events with too much activity in the opposite edge region
  if ( (nr->eventq.oppEdgeStrips>3)
       && (nr->eventq.oppEdgePH>1000)
       && (nr->srevent.energyGeV<5.0 )
         && (nr->srshower.planes>nr->srtrack.planes) ) // don't cut out CC
    return false;
                                                                                
  // make additional deltaZ cuts if (40ns<|minDeltaT|<120ns)
  if (TMath::Abs(nr->eventq.closeTimeDeltaZ)<1.0
      && TMath::Abs(nr->eventq.minTimeSeparation)<120e-9)
    return false;
                                                                                
  // Survived all cuts
  return true;
}

void NueStandard::SetDefaultOscParam()
{
  Double_t dm2_12 = 8.0e-5;  // best fit SNO
  Double_t dm2_23 = 2.43e-3;

  Double_t par[9] = {0};
  par[OscPar::kL] = 735.0;
  par[OscPar::kTh23] = 3.1415926/4.0;
  par[OscPar::kTh12] = 0.59365;    // Sin2(2Th12) = 0.86
  par[OscPar::kTh13] = 0.19885;  // Sin2(2Th13) = 0.15;
  par[OscPar::kDeltaM23] = dm2_23; //normal heirarchy
  par[OscPar::kDeltaM12] = dm2_12;
  par[OscPar::kDensity] = 2.75; //standard rock density
  par[OscPar::kDelta] = 0;
  par[OscPar::kNuAntiNu] = 1;

  fOscGen.SetOscParam(par);
}

void NueStandard::SetDefaultOscParamNoNue()
{
  NueStandard::SetDefaultOscParam();
  NueStandard::SetOscParam(OscPar::kTh13, 0);
}

void NueStandard::SetOscParamBestFitANN()
{
  NueStandard::SetDefaultOscParam();
  NueStandard::SetOscParam(OscPar::kTh13, 0.172088);
}

void NueStandard::SetOscParamBestFitkNN()
{ 
  NueStandard::SetDefaultOscParam(); 
  NueStandard::SetOscParam(OscPar::kTh13, 0.145518);
} 

void NueStandard::SetOscNoMatter()
{
   fOscGen.SetOscParam(OscPar::kDensity, 1e-9);
}

void NueStandard::SetOscParam(double *par)
{
   fOscGen.SetOscParam(par);
}

void NueStandard::SetOscParam(OscPar::OscPar_t pos, double val)
{
  fOscGen.SetOscParam(pos, val);
}

double NueStandard::GetOscWeight(int nuFlavor, int nonOsc, double E)
{
   return fOscGen.Oscillate(nuFlavor, nonOsc, E);
}

void NueStandard::FillDefaultOscParam(double* par)
{
  Double_t dm2_12 = 8.0e-5;  // best fit SNO
  Double_t dm2_23 = 2.43e-3;
                                                                                      
  par[OscPar::kL] = 735.0;
  par[OscPar::kTh23] = 3.1415926/4.0;
  par[OscPar::kTh12] = 0.59365;    // Sin2(2Th12) = 0.86
  par[OscPar::kTh13] = 0.19885;  // Sin2(2Th13) = 0.15;
  par[OscPar::kDeltaM23] = dm2_23; //normal heirarchy
  par[OscPar::kDeltaM12] = dm2_12;
  par[OscPar::kDensity] = 2.75; //standard rock density
  par[OscPar::kDelta] = 0;
  par[OscPar::kNuAntiNu] = 1;
}

void NueStandard::GetOscParam(double* par){
  fOscGen.GetOscParam(par);
}

Bool_t NueStandard::IsRun1(NueRecord* nr)
{
 return(nr->GetHeader().GetVldContext().GetTimeStamp().GetSec() < 1145036420);
}


Bool_t NueStandard::IsRun2(NueRecord* nr)
{
 return(nr->GetHeader().GetVldContext().GetTimeStamp().GetSec() > 1145036420 && nr->GetHeader().GetVldContext().GetTimeStamp().GetSec() < 1190028111);
}

Bool_t NueStandard::IsRun3(NueRecord* nr)
{
 return(nr->GetHeader().GetVldContext().GetTimeStamp().GetSec() > 1190028111 && nr->GetHeader().GetVldContext().GetTimeStamp().GetSec() < 1249594050);
}


 Double_t NueStandard::GetMCWeights(NueRecord *nr)
 {
  //Function to calculate MC event weights to correct for various effect
   Double_t eventWeight = 1.0;
   
  //Check if it is MC
   if(nr->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)
   {
    //Correct MC flux (i.e. SKZP and Helium) 
    eventWeight *= NueStandard::GetSKZPBeamWeight(nr);
   }

   return(eventWeight);   
 }



 Double_t NueStandard::GetSKZPBeamWeight(NueRecord *nr)
 {
  //One time intialization
   static SKZPWeightCalculator skzpWC("DetXs",true);  //static because we only want this created once
                                                       //and not every time it's called within a loop
   static bool firstTimeFunctionCalled = true;
  
   if(firstTimeFunctionCalled)
   {
    //print out the configuration information once
     skzpWC.PrintReweightConfig(std::cout);
     firstTimeFunctionCalled = false;
   }

  //Now get the weight
   Double_t eventWeight = 1.0;
   
  //Check if it is MC
   if(nr->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)
   {

    //Correct MC flux (i.e. SKZP and Helium) 
    double pt = (nr->fluxinfo.tpx * nr->fluxinfo.tpx);
           pt += (nr->fluxinfo.tpy * nr->fluxinfo.tpy);
           pt = sqrt(pt);

    eventWeight = skzpWC.GetBeamWeight(
                                nr->GetHeader().GetVldContext().GetDetector() /*int det*/,
                                BeamType::ToZarko( nr->GetHeader().GetBeamType() ) /*int Ibeam*/,

                                nr->fluxinfo.tptype /*int tptype*/,
                                pt /*double pt*/,
                                nr->fluxinfo.tpz /*double pz*/,
                      
                                nr->mctrue.nuEnergy /*double true_enu*/,
                                nr->mctrue.nonOscNuFlavor /*int inu*/,  //want weight before neutrino oscillations
                                nr->GetHeader().GetVldContext() /*VldContext vc*/
                               );
   }

   return(eventWeight);   
 }

void NueStandard::ModifyANNPID(NueRecord *nr) {
 
 // This is only a temporary solution for the ANN14 PID.
 
  static bool firstcall = true;
  static TMultiLayerPerceptron* fneuralNet_14inp_daikon04 = 0;

  double ann14pid = -9999.9;

  char *srt_dir = getenv("SRT_PRIVATE_CONTEXT");
  char annfile[1000];

  if (firstcall) {
    srt_dir = getenv("SRT_PRIVATE_CONTEXT");
    sprintf(annfile,"%s/NueAna/data/ann14_400_14_9.root",srt_dir);
    ifstream Test2(annfile);
    if (!Test2){
      srt_dir = getenv("SRT_PUBLIC_CONTEXT");
      sprintf(annfile,"%s/NueAna/data/ann14_400_14_9.root",srt_dir);
      ifstream Test_again2(annfile);
      if (!Test_again2){
        cout<<"Couldn't find ANN object, blame Jiajie"<<endl;
        exit(0);
      }
    }

    static TFile *f = TFile::Open(annfile);
    fneuralNet_14inp_daikon04 = (TMultiLayerPerceptron*) f->Get("mlp");
  }

  if (!fneuralNet_14inp_daikon04) {
    cout << "Couldn't find ANN object, blame Jiajie" << endl;
    exit(0);
  }

  if (firstcall) {
    cout << "Get ann14 value from : " << annfile << endl;
    firstcall = false;
  }

  NueConvention::NueEnergyCorrection(nr); 

  double paras[14];
  paras[0] = nr->shwfit.par_a;
  paras[1] = nr->shwfit.par_b;
  paras[2] = nr->shwfit.uv_rms_9s_2pe_dw;
  paras[3] = nr->fracvars.fract_road;
  paras[4] = nr->shwfit.LongE;
  paras[5] = nr->fracvars.shw_max;
  paras[6] = nr->fracvars.shw_slp;
  paras[7] = nr->srevent.phNueGeV;
  paras[8] = nr->fracvars.dis2stp;
  paras[9] = nr->fracvars.fract_1_plane;
  paras[10] = nr->fracvars.fract_5_planes;
  paras[11] = nr->fracvars.fract_6_counters;
  paras[12] = nr->fracvars.fract_20_counters;
  paras[13] = TMath::Abs(nr->srevent.endPlane - nr->srevent.begPlane);

  int pass = 1;

  if (nr->shwfit.par_a<-1000) pass = 0;
  if (nr->shwfit.par_b<-1000) pass = 0;
  if (nr->shwfit.uv_rms_9s_2pe_dw<-1000) pass = 0;
  if (nr->fracvars.fract_1_plane<-1000) pass = 0;
  if (nr->fracvars.fract_5_planes<-1000) pass = 0;
  if (nr->fracvars.fract_6_counters<-1000) pass = 0;
  if (nr->fracvars.fract_20_counters<-1000) pass = 0;
  if (nr->fracvars.fract_road<-1000) pass = 0;
  if (nr->shwfit.LongE<-1000) pass = 0;
  if (nr->shwfit.LongE>1000) pass = 0;
  if (nr->fracvars.shw_max<-1000) pass = 0;
  if (nr->fracvars.dis2stp<-1000) pass = 0;
  if (nr->fracvars.shw_slp<0) pass = 0;
  if (nr->srevent.endPlane<-1000) pass = 0;
  if (nr->srevent.begPlane<-1000) pass = 0;
  
  if (pass) {
    ann14pid = fneuralNet_14inp_daikon04->Evaluate(0, paras);
  } 

  // set the ann.pid value to be ANN14 
  nr->ann.pid = ann14pid;
}

---