MadTVAnalysis.cxx

Go to the documentation of this file.

//
// MadTVAnalysis
//
// Tricia's MadAnalysis.  Copied from the Feb 3 version of the jan06 branch of 
// MadMKAnalysis
//

#include <set>
#include <vector>
#include "TSystem.h"
#include "TClonesArray.h"
#include "Mad/MadTVAnalysis.h"
#include "Mad/NearbyEvents.h"
#include "Mad/Moments.h"
#include "DataUtil/EnergyCorrections.h"
// tricia's beam stuff
//#include "Mad/BeamMonMap.h"
#include "MCReweight/NuParent.h"
#include "MCReweight/BeamEnergyCalculator.h"
//#include "MCReweight/BeamType.h"
#include "Conventions/BeamType.h"
#include "Conventions/ReleaseType.h"
#include "DcsUser/CoilTools.h"
#include "DataUtil/DataQualDB.h"

#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "SpillTiming/SpillTimeFinder.h"

#include "DatabaseInterface/DbiResultPtr.h"
#include "DcsUser/Dcs_Mag_Near.h"

#include "PhysicsNtuple/Store/Interface.h"
#include "NtupleUtils/LISieve.h"

// david's list of bad run times
#include "Mad/fddataquality.h"

const Float_t MadTVAnalysis::kVetoEndZ = 1.2154;
const Float_t MadTVAnalysis::kTargEndZ = 3.5914;
const Float_t MadTVAnalysis::kCalorEndZ = 7.1554;
const Float_t MadTVAnalysis::kHalfCalorEndZ = kTargEndZ+(kCalorEndZ-kTargEndZ)*0.5;
const Float_t MadTVAnalysis::kSpecEndZ = 16.656;

//June 2007, beam spot update after survey
//const Float_t MadTVAnalysis::kXcenter = 1.4885;
//const Float_t MadTVAnalysis::kYcenter = 0.1397;
const Float_t MadTVAnalysis::kXcenter = 1.4828;
const Float_t MadTVAnalysis::kYcenter = 0.2384;


const Double_t MadTVAnalysis::fEarlyActivityWindowSize=10.;
const Double_t MadTVAnalysis::fPMTTimeConstant=700.;
const Int_t MadTVAnalysis::fNPlaneEarlyAct=30;
const Double_t MadTVAnalysis::fEarlyPlaneSphere=0.5;

MadTVAnalysis::MadTVAnalysis(TChain *chainSR,TChain *chainMC,
                             TChain *chainTH,TChain *chainEM)
  : MadAnalysis(chainSR, chainMC, chainTH, chainEM)
{
  openedabpidfile=false;
  fMCBeam=BeamType::kUnknown;
  fDataUseMCBeam=false;
  if(!chainSR) {
    record = 0;
    emrecord = 0;
    mcrecord = 0;
    threcord = 0;
    Clear();
    whichSource = -1;
    isMC = true;
    isTH = true;
    isEM = true;
    Nentries = -1;
    return;
  }
  static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
  fBECConfig.Set("beam:flux_file",default_bec_file);
  fBECConfig.Set("beam:beam_type",BeamType::kLE);
  fBECConfig.Set("beam:detector",Detector::kNear);
  fBECConfig.Set("beam:low_energy_limit",0.5);
  fBECConfig.Set("beam:high_energy_limit",60.0);

  
  //  InitChain(chainSR,chainMC,chainTH,chainEM);
  whichSource = 0;  

  phnti = new Anp::Interface();

}

//********************************************************
MadTVAnalysis::MadTVAnalysis(JobC *j,string path,int entries) : MadAnalysis(j, path, entries)
{
  openedabpidfile=false;
  fMCBeam=BeamType::kUnknown;
  fDataUseMCBeam=false;

  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  Clear();
  isMC = true;
  isTH = true;
  isEM = true;
  Nentries = entries;
  whichSource = 1;
  jcPath = path;
  fJC = j;
  fChain = NULL;

  
  static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
  fBECConfig.Set("beam:flux_file",default_bec_file);
  fBECConfig.Set("beam:beam_type",BeamType::kLE);
  fBECConfig.Set("beam:detector",Detector::kNear);
  fBECConfig.Set("beam:low_energy_limit",0.5);
  fBECConfig.Set("beam:high_energy_limit",60.0);

  phnti = new Anp::Interface();

}

//********************************************************
MadTVAnalysis::~MadTVAnalysis()
{
  if(phnti){
    delete phnti;
    phnti=0;
  }
}

void MadTVAnalysis::CreatePAN(std::string tag, const char* /*bmonpath*/)
{
  
  // is this MC??
  this->GetEntry(0);
  const bool file_is_mc
    =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);

  //set mc and recoversions for pid histogram retrevial
  //  string recover = EnergyCorrections::GetCorrectionAsString();
  string mcver = "daikon"; //note this is hard coded and bad, but I dont' care.
  ReleaseType::Release_t reltype = ReleaseType::MakeReleaseType(strecord->GetTitle(),mcver);
  string recover = ReleaseType::AsString(reltype);
  
  cout<<"Running pans over "<<recover<<" "<<mcver<<endl;

  //updated ireg to from rustems for his cc-pid selection 12-21-2007
  //set up anti numu selection and rustem's cc pid
  Registry ireg(false);
  
  ireg.Set("InterfaceConfigPath", rustemconfigfile.c_str());
  if(ntpHeader->GetVldContext().GetDetector() == Detector::kNear)
  {
    ireg.Set("FillkNNFilePath", rustempidfile_near.c_str());
  }
  else if(ntpHeader->GetVldContext().GetDetector() == Detector::kFar)
  {
    ireg.Set("FillkNNFilePath", rustempidfile_far.c_str());
  }
  
  //
  // Do not erase events that fail SelectAntiNeutrino algorithm
  //
  ireg.Set("QuietInterface", "yes");
  ireg.Set("SelectAntiNeutrinoErase", "no");
  
  ireg.LockValues();
  phnti->Config(ireg);

  //PAN Quantities 
  //Truth:
  Float_t true_enu;       // true neutrino energy (GeV)
  Float_t true_emu;       // true muon energy
  Float_t true_eshw;      // true shower energy
  Float_t true_x;         // true x
  Float_t true_y;         // true y
  Float_t true_q2;        // true q2
  Float_t true_w2;        // true w2
  Float_t true_dircosneu; // true muon dircos w.r.t neutrino
  Float_t true_dircosz;   // track z-direction cosine
  Float_t true_vtxx;      // true x vtx
  Float_t true_vtxy;      // true y vtx
  Float_t true_vtxz;      // true z vtx
  Int_t true_pitt_fid; // was event truly in fiducial volue?
  Float_t true_trk_cmplt;  //true track completeness

  //For Neugen:
  Float_t flavour;        // true flavour: 1-e 2-mu 3-tau
  Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
  Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
  Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
  Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
  Int_t inu;              // stdhep inu
  Int_t inunoosc;              // stdhep inu unoscilllated
  Int_t resnum; // IdHEP%1000 of resonance state

  Short_t npp; // # final state pi+
  Short_t npm; // # fs pi-
  Short_t npz; // # fs pi0
  Short_t np; // # fs p
  Short_t nn; // # fs n
  float epp; // energy final state pi+
  float epm; // energy fs pi-
  float epz; // energy fs pi0
  float ep; // energy fs p
  float en; // energy fs n

  
  //Reco Quantities
  Int_t detector;         // Near = 1, Far = 2;
  Int_t run;              // run number
  Int_t snarl;            // snarl number
  Int_t subrun;           // subrun number
  UInt_t trgsrc;          // trigger source
  double trgtime;         // trigger time
  Double_t vld_sec;          // seconds field of header vld context
  Int_t spilltype;        // comes from mdSpillTypeEnum
                          // 0=none, 1=reported, 2=predicted
                          // 3=fake, 4=local

  Int_t nevent;           // number of events in snarl
  Int_t event;            // event index
  Int_t mcevent;          // mc event index
  Int_t ntrack;           // number of tracks in event
  Int_t nshower;          // number of showers in event
  Int_t nstrip;           // number of strips in event

  Int_t is_fid;           // pass fid cut. true = 1, false = 0
  Int_t is_fid_2007;      // pass fid cut. true = 1, false = 0
  Int_t is_cev;           // fully contained. true = 1, false = 0 
  Int_t sr_contained;     // containment support from reconstruction. true=1,false=0 
  Int_t is_mc;            // is a corresponding mc event found
  Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
  Int_t is_pitt_fid;      // passes pitt ND vtx fid cut true=1, false=0
  Int_t is_trk_fid;       // is the trkvtx in fid (pitt fid for ND)
                          // true=1, false=0, notrack = -1
  Int_t nd_radial_fid;    // a bitfield, bits correspond to r cuts
  Int_t pitt_evt_class;   // CC event class as defined by pitt group
  // 1 = track is fully contained in the upstream region
  // 2 = track is partially contained in the upstream region
  // 3 = track is fully contained in the downstream region
  // 4 = track is partially contained in the downstream region
  Int_t pitt_trk_qual;   // pitt tracking quality cuts
  Int_t duvvtx;          //delta (Uvtx-Vvtx)

  Float_t pid0;           // pid parameter (usual method)
  Float_t pid1;           // pid parameter (alternative method 1)
  Float_t pid2;           // pid parameter (alternative method 2)
  Int_t pass_track;       // the primary track passes track cuts
  Int_t trk_fit_pass;       // trk.fit.pass
  Int_t pass;             // pass analysis cuts. true = 1, false = 0

  Int_t emu_meth;         // method used to determine muon energy
  Float_t reco_enu;       // reco neutrino energy
  Float_t reco_emu;       // reco muon energy
  Float_t reco_mushw;     // showers along muon track
  Float_t reco_eshw;      // reco shower energy (largest shower
  Float_t reco_wtd_eshw;    // as above but weighted 
  Float_t reco_vtx_eshw;    // reco shower energy (=0 if no vertex shower)
  Float_t reco_vtx_wtd_eshw;// as above but weighted (=0 if no vertex shower)
  
  Float_t reco_eshw_uncorr; // uncorrected shower energy
  Float_t reco_wtd_eshw_uncorr; // uncorrected weighted shower energy

  Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
  Float_t reco_dircosx;   // reco track vtx x-dircos
  Float_t reco_dircosy;   // reco track vtx y-dircos
  Float_t reco_dircosz;   // reco track vtx z-dircos

  Float_t reco_vtxx;      // reco x vtx
  Float_t reco_vtxy;      // reco y vtx
  Float_t reco_vtxz;      // reco z vtx
  Float_t reco_vtxr;// radial position of vertex

  Float_t reco_trk_vtxx;      // reco x track vtx
  Float_t reco_trk_vtxy;      // reco y track vtx
  Float_t reco_trk_vtxz;      // reco z track vtx
  Float_t reco_trk_vtxr;// radial position of vertex

  Float_t reco_shw_vtxx;      // reco x track vtx
  Float_t reco_shw_vtxy;      // reco y track vtx
  Float_t reco_shw_vtxz;      // reco z track vtx
  Float_t reco_shw_vtxr;// radial position of vertex

  Float_t reco_trk_endx;      // reco x track end
  Float_t reco_trk_endy;      // reco y track end
  Float_t reco_trk_endz;      // reco z track end
  Float_t reco_trk_endr;// radial position of vertex

  Int_t evt_nplanes;        //number of planes in event above some ph
  Int_t trk_nplanes;        //number of planes in track above some ph
  Int_t slc_nplanes;        //number of planes in slice above some ph
  Int_t shw_nplanes_cut;        //number of planes in shower above some ph
  Int_t shw_nplanes_raw;        //number of planes in shower
  Float_t shw_ph_cut, shw_ph_raw; // sigcor in shower, w/ & w/o ph cut
  
  Int_t evt_nstrips;        //number of strips in event above some ph
  Int_t trk_nstrips;        //number of strips in track above some ph
  Int_t slc_nstrips;        //number of strips in slice above some ph


  // reconstructed kinematics
  // formulae given for high-energy CC interactions
  // Assume nu = Ehad
  Float_t reco_y;// y = (Enu-Emu)/Enu
  Float_t reco_Q2;// Q2 = -q^{2} = 2 Enu Emu (1 - cos(theta_mu))
  Float_t reco_x;// x = Q2 / (2M* Ehad)  {M=nucleon mass}
  Float_t reco_W2;// W2 = M^{2} + q^{2} + 2M(Enu-Emu)

  // reco quantities, assuming a QE event
  Float_t reco_qe_enu;    // reco qe neutrino energy
  Float_t reco_qe_q2;     // reco qe q2

  
  Float_t evtlength;      // reco event length
  Float_t shwlength;      // reco shower length
  Float_t trklength;      // reco track length (planes)
  Int_t ntrklike;      // number of track-like planes
  Float_t trkmomrange;    // reco track momentum from range
  Float_t trkqp;          // reco track q/p
  Float_t trkeqp;         // reco track q/p error
  Float_t trkphfrac;      // reco pulse height fraction in track
  Float_t trkphplane;     // reco average track pulse height/plane
  Float_t trkchi2;
  Int_t trkndf;


  Int_t trkend, trkendu, trkendv; // track end plane, u, v
  Int_t shwend,shwbeg; // shower begin,end planes

  Int_t nss[6];//
  Float_t ess[6];//
  Int_t ntotss;



  // variables to tag nearest event  in space and time
  // nvsi = nearest vertex space : index
  // nvsd = nearest vertex space : distance
  // nvst = nearest vertex space : distance
  // nvsp = nearest vertex space : pulse height
  // nvsevt = nearest vertex space : mc event
  // nvti = nearest vertex time : index
  // ....  
  // bvsi = back vertex space : index
  // the "back" variables hold the same info as the "near" variables
  // but for that nearest event
  // I'm interested in this case
  //  
  //  (this vtx)------>(nearest)--->(nearest to it)
  //
  // if you are going to try to recombine events
  // you need to somehow account for this
  //

  Int_t nvsi, bvsi, nvti, bvti;
  Int_t nvsevt, bvsevt, nvtevt, bvtevt;
  Float_t nvsd, bvsd, nvtd, bvtd;
  Float_t nvst, bvst, nvtt, bvtt;
  Float_t nvsp, bvsp, nvtp, bvtp;

  // this event's pulseheight (in sigcor)
  Float_t evtph;
  Float_t evtphgev;

  // the pulse height of the event and track
  // in the fully and partially instrumented regions
  Float_t evtsigfull,trksigfull,evtsigpart,trksigpart;


  //time structure of event
  double evttimemin, evttimemax, evttimeln;
  double evttimeminphc, evttimemaxphc, evttimelnphc;

  //early activity variables
  float earlywph, earlywphpw, earlywphsphere;
  float ph1000, ph500, ph100;
  float lateph1000, lateph500, lateph100, lateph50;

  //duplicate reconstructed strip variables
  int nduprs; //number of reco strips that are duplicated in another event
  float dupphrs; //pulseheight of reco strips 
                 //that are duplicated in another event

  //time last li event in snarl, -1 if none;
  double litime;
  int is_li_sieve=0;

  //demux failures (for fd)
  UChar_t dmx_stat;
  //  std::cout<<"Second Print"<<std::endl;

  // david's FD data quality
  Int_t dp_fd_quality;

  //brian's cut on lowph ratio to remove junk events 
  //in FD that pass the dmx_stat variable
  float lowphrat;

  //Trees
  // output file
  std::string savename = "PAN_" + tag + ".root";
  TFile save(savename.c_str(),"RECREATE"); 
  save.cd();
  TTree *tree = new TTree("pan","pan");
  //Truth
  tree->Branch("true_enu",&true_enu,"true_enu/F",32000);
  tree->Branch("true_emu",&true_emu,"true_emu/F",32000);
  tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000);
  tree->Branch("true_x",&true_x,"true_x/F",32000);
  tree->Branch("true_y",&true_y,"true_y/F",32000);
  tree->Branch("true_q2",&true_q2,"true_q2/F",32000);
  tree->Branch("true_w2",&true_w2,"true_w2/F",32000);
  Float_t nu_px, nu_py, nu_pz;
  Float_t tar_px, tar_py, tar_pz, tar_e;
  tree->Branch("nu_px", &nu_px, "nu_px/F");
  tree->Branch("nu_py", &nu_py, "nu_py/F");
  tree->Branch("nu_pz", &nu_pz, "nu_pz/F");
  tree->Branch("tar_px", &tar_px, "tar_px/F");
  tree->Branch("tar_py", &tar_py, "tar_py/F");
  tree->Branch("tar_pz", &tar_pz, "tar_pz/F");
  tree->Branch("tar_e", &tar_e, "tar_e/F");


  tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000);
  tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000);
  tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000);
  tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000);
  tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000);
  tree->Branch("true_pitt_fid",&true_pitt_fid,"true_pitt_fid/I",32000);
  tree->Branch("true_trk_cmplt",&true_trk_cmplt,"true_trk_cmplt/F",32000);

  //Neugen
  tree->Branch("flavour",&flavour,"flavour/F",32000);
  tree->Branch("process",&process,"process/I",32000);
  tree->Branch("nucleus",&nucleus,"nucleus/I",32000);
  tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000);
  tree->Branch("inu",&inu,"inu/I",32000);
  tree->Branch("inunoosc",&inunoosc,"inunoosc/I",32000);

  tree->Branch("initial_state",&initial_state,"initial_state/I",32000);
  tree->Branch("resnum",&resnum,"resnum/I",32000);
  tree->Branch("npp",&npp,"npp/S",32000);
  tree->Branch("npm",&npm,"npm/S",32000);
  tree->Branch("npz",&npz,"npz/S",32000);
  tree->Branch("np",&np,"np/S",32000);
  tree->Branch("nn",&nn,"nn/S",32000);

  tree->Branch("epp",&epp,"epp/F",32000);
  tree->Branch("epm",&epm,"epm/F",32000);
  tree->Branch("epz",&epz,"epz/F",32000);
  tree->Branch("ep",&ep,"ep/F",32000);
  tree->Branch("en",&en,"en/F",32000);

  //Reco
  tree->Branch("detector",&detector,"detector/I",32000);
  tree->Branch("run",&run,"run/I",32000);
  tree->Branch("snarl",&snarl,"snarl/I",32000);
  tree->Branch("subrun",&subrun,"subrun/I",32000);
  tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000);
  tree->Branch("trgtime",&trgtime,"trgtime/D",32000);
  tree->Branch("vld_sec",&vld_sec,"vld_sec/D",32000);


  tree->Branch("spilltype",&spilltype,"spiltype/I",32000);
  tree->Branch("nevent",&nevent,"nevent/I",32000);
  tree->Branch("event",&event,"event/I",32000);
  tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
  tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
  tree->Branch("nshower",&nshower,"nshower/I",32000);
  tree->Branch("nstrip",&nstrip,"nstrip/I",32000);
  
  tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
  tree->Branch("is_fid_2007",&is_fid_2007,"is_fid_2007/I",32000);
  tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
  tree->Branch("sr_contained",&sr_contained,"sr_contained/I",32000);
  tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
  tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
  tree->Branch("pass_track",&pass_track,"pass_track/I",32000);
  tree->Branch("trk_fit_pass",&trk_fit_pass,"trk_fit_pass/I",32000);
  tree->Branch("emu_meth",&emu_meth,"emu_meth/I",32000);
  tree->Branch("is_pitt_fid",&is_pitt_fid,"is_pitt_fid/I",32000);
  tree->Branch("is_trk_fid",&is_trk_fid,"is_trk_fid/I",32000);

  tree->Branch("nd_radial_fid",&nd_radial_fid,"nd_radial_fid/I",32000);
  tree->Branch("pitt_evt_class",&pitt_evt_class,"pitt_evt_class/I",32000);
  tree->Branch("pitt_trk_qual",&pitt_trk_qual,"pitt_trk_qual/I",32000);
  tree->Branch("duvvtx",&duvvtx,"duvvtx/I",32000);

  tree->Branch("pid0",&pid0,"pid0/F",32000);
  tree->Branch("pid1",&pid1,"pid1/F",32000);
  tree->Branch("pid2",&pid2,"pid2/F",32000);
  tree->Branch("pass",&pass,"pass/I",32000);

  float med_qe_pid;
  int med_qe_pass;
  tree->Branch("med_qe_pid",&med_qe_pid,"med_qe_pid/F",32000);
  tree->Branch("med_qe_pass",&med_qe_pass,"med_qe_pass/I",32000);

  double niki_cc_pid;
  tree->Branch("niki_cc_pid",&niki_cc_pid,"niki_cc_pid/D",32000);
  float dave_cc_pid;
  tree->Branch("dave_cc_pid",&dave_cc_pid,"dave_cc_pid/F",32000);


  tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
  tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
  tree->Branch("reco_mushw",&reco_mushw,"reco_mushw/F",32000);
  tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
  tree->Branch("reco_wtd_eshw",&reco_wtd_eshw,"reco_wtd_eshw/F",32000);
  tree->Branch("reco_vtx_eshw",&reco_vtx_eshw,"reco_vtx_eshw/F",32000);
  tree->Branch("reco_vtx_wtd_eshw",&reco_vtx_wtd_eshw,"reco_vtx_wtd_eshw/F",32000);
  tree->Branch("reco_eshw_uncorr",&reco_eshw_uncorr,"reco_eshw_uncorr/F",32000);
  tree->Branch("reco_wtd_eshw_uncorr",&reco_wtd_eshw_uncorr,"reco_wtd_eshw_uncorr/F",32000);
 
  tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
  tree->Branch("reco_dircosx",&reco_dircosx,"reco_dircosx/F",32000);
  tree->Branch("reco_dircosy",&reco_dircosy,"reco_dircosy/F",32000);
  tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);

  tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
  tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
  tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
  tree->Branch("reco_vtxr",&reco_vtxr,"reco_vtxr/F",32000);

  tree->Branch("reco_trk_vtxx",&reco_trk_vtxx,"reco_trk_vtxx/F",32000);
  tree->Branch("reco_trk_vtxy",&reco_trk_vtxy,"reco_trk_vtxy/F",32000);
  tree->Branch("reco_trk_vtxz",&reco_trk_vtxz,"reco_trk_vtxz/F",32000);
  tree->Branch("reco_trk_vtxr",&reco_trk_vtxr,"reco_trk_vtxr/F",32000);

  tree->Branch("reco_shw_vtxx",&reco_shw_vtxx,"reco_shw_vtxx/F",32000);
  tree->Branch("reco_shw_vtxy",&reco_shw_vtxy,"reco_shw_vtxy/F",32000);
  tree->Branch("reco_shw_vtxz",&reco_shw_vtxz,"reco_shw_vtxz/F",32000);
  //  tree->Branch("reco_shw_vtxr",&reco_shw_vtxr,"reco_shw_vtxr/F",32000);

  tree->Branch("reco_trk_endx",&reco_trk_endx,"reco_trk_endx/F",32000);
  tree->Branch("reco_trk_endy",&reco_trk_endy,"reco_trk_endy/F",32000);
  tree->Branch("reco_trk_endz",&reco_trk_endz,"reco_trk_endz/F",32000);
  tree->Branch("reco_trk_endr",&reco_trk_endr,"reco_trk_endr/F",32000);

  tree->Branch("evt_nplanes",&evt_nplanes,"evt_nplanes/I",32000);
  tree->Branch("trk_nplanes",&trk_nplanes,"trk_nplanes/I",32000);
  tree->Branch("slc_nplanes",&slc_nplanes,"slc_nplanes/I",32000);
  tree->Branch("shw_nplanes_cut",&shw_nplanes_cut,"shw_nplanes_cut/I",32000);
  tree->Branch("shw_nplanes_raw",&shw_nplanes_raw,"shw_nplanes_raw/I",32000);

  tree->Branch("shw_ph_cut",&shw_ph_cut,"shw_ph_cut/F",32000);
  tree->Branch("shw_ph_raw",&shw_ph_raw,"shw_ph_raw/F",32000);


  tree->Branch("evt_nstrips",&evt_nstrips,"evt_nstrips/I",32000);
  tree->Branch("trk_nstrips",&trk_nstrips,"trk_nstrips/I",32000);
  tree->Branch("slc_nstrips",&slc_nstrips,"slc_nstrips/I",32000);

  tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
  tree->Branch("reco_Q2",&reco_Q2,"reco_Q2/F",32000);
  tree->Branch("reco_x",&reco_x,"reco_x/F",32000);
  tree->Branch("reco_W2",&reco_W2,"reco_W2/F",32000);

  tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
  tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);

  tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
  tree->Branch("trklength",&trklength,"trklength/F",32000);
  tree->Branch("shwlength",&shwlength,"shwlength/F",32000);
  tree->Branch("ntrklike",&ntrklike,"ntrklike/I",32000);
  tree->Branch("trkend",&trkend,"trkend/I",32000);
  tree->Branch("trkendu",&trkendu,"trkendu/I",32000);
  tree->Branch("trkendv",&trkendv,"trkendv/I",32000);
  tree->Branch("shwbeg",&shwbeg,"shwbeg/I",32000);
  tree->Branch("shwend",&shwend,"shwend/I",32000);


  Int_t trknstp;
  tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
  tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
  tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
  tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
  tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
  tree->Branch("trkchi2",&trkchi2,"trkchi2/F",32000);
  tree->Branch("trkndf",&trkndf,"trkndf/I",32000);
  tree->Branch("trknstp",&trknstp,"trknstp/I",32000);


  // vertex track back variables (described above and in code)
  tree->Branch("nvsi", &nvsi, "nvsi/I");
  tree->Branch("nvsevt", &nvsevt, "nvsevt/I");
  tree->Branch("nvsd", &nvsd, "nvsd/F");
  tree->Branch("nvst", &nvst, "nvst/F");
  tree->Branch("nvsp", &nvsp, "nvsp/F");
  tree->Branch("nvti", &nvti, "nvti/I");
  tree->Branch("nvtevt", &nvtevt, "nvtevt/I");
  tree->Branch("nvtd", &nvtd, "nvtd/F");
  tree->Branch("nvtt", &nvtt, "nvtt/F");
  tree->Branch("nvtp", &nvtp, "nvtp/F");
  /*
  tree->Branch("bvsi", &bvsi, "bvsi/I");
  tree->Branch("bvsevt", &bvsevt, "bvsevt/I");
  tree->Branch("bvsd", &bvsd, "bvsd/F");
  tree->Branch("bvst", &bvst, "bvst/F");
  tree->Branch("bvsp", &bvsp, "bvsp/F");
  tree->Branch("bvti", &bvti, "bvti/I");
  tree->Branch("bvtevt", &bvtevt, "bvtevt/I");
  tree->Branch("bvtd", &bvtd, "bvtd/F");
  tree->Branch("bvtt", &bvtt, "bvtt/F");
  tree->Branch("bvtp", &bvtp, "bvtp/F");
  */

  tree->Branch("evtph", &evtph, "evtph/F");
  tree->Branch("evtphgev", &evtphgev, "evtphgev/F");
  
  tree->Branch("evtsigfull", &evtsigfull, "evtsigfull/F");
  tree->Branch("evtsigpart", &evtsigpart, "evtsigpart/F");

  tree->Branch("trksigfull", &trksigfull, "trksigfull/F");
  tree->Branch("trksigpart", &trksigpart, "trksigpart/F");

  tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
  tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
  tree->Branch("evttimeln",&evttimeln,"evttimeln/D");

  tree->Branch("evttimeminphc",&evttimeminphc,"evttimeminphc/D");
  tree->Branch("evttimemaxphc",&evttimemaxphc,"evttimemaxphc/D");
  tree->Branch("evttimelnphc",&evttimelnphc,"evttimelnphc/D");

  tree->Branch("earlywph",&earlywph,"earlywph/F");
  tree->Branch("earlywphpw",&earlywphpw,"earlywphpw/F");
  tree->Branch("earlywphsphere",&earlywphsphere,"earlywphsphere/F");
  tree->Branch("ph1000",&ph1000,"ph1000/F");
  tree->Branch("ph500",&ph500,"ph500/F");
  tree->Branch("ph100",&ph100,"ph100/F");
  tree->Branch("lateph1000",&lateph1000,"lateph1000/F");
  tree->Branch("lateph500",&lateph500,"lateph500/F");
  tree->Branch("lateph100",&lateph100,"lateph100/F");
  tree->Branch("lateph50",&lateph50,"lateph50/F");
  Int_t largest_event;
  tree->Branch("largest_event",&largest_event,"largest_event/I");

  tree->Branch("nduprs",&nduprs,"ndubrs/I");
  tree->Branch("dupphrs",&dupphrs,"dupphrs/F");

  tree->Branch("nss", nss, "nss[6]/I");
  tree->Branch("ess", ess, "ess[6]/F");
  tree->Branch("ntotss", &ntotss, "ntotss/I");

  Float_t etu, etv, elu, elv;
  tree->Branch("etu", &etu, "etu/F"); // transveres u and v energy
  tree->Branch("etv", &etv, "etv/F"); // uses subshowers
  tree->Branch("elu", &elu, "elu/F"); // longitudinal energy
  tree->Branch("elv", &elv, "elv/F");

  Float_t swu,swv,wswu,wswv; 
  tree->Branch("swu",&swu, "swu/F");// shower width in u
  tree->Branch("swv",&swv, "swv/F");// shower width in u
  tree->Branch("wswu",&wswu, "wswu/F");// shower width in u
  tree->Branch("wswv",&wswv, "wswv/F");// shower width in u

  // beam energy weights
  // 1/E flux, le, z=050, z=100, z=200, z=250 cm
  Float_t bec_inve, bec_le, bec_050, bec_100, bec_200, bec_250;
  
  tree->Branch("bec_inve", &bec_inve, "bec_inve/F");
  tree->Branch("bec_le", &bec_le, "bec_le/F");
  tree->Branch("bec_050", &bec_050, "bec_050/F");
  tree->Branch("bec_100", &bec_100, "bec_100/F");
  tree->Branch("bec_200", &bec_200, "bec_200/F");
  tree->Branch("bec_250", &bec_250, "bec_250/F");

  //li time
  tree->Branch("litime",&litime,"litime/D");
  Int_t rc_boundary;
  tree->Branch("rc_boundary",&rc_boundary,"rc_boundary/I");
  tree->Branch("is_li_sieve",&is_li_sieve,"is_li_sieve/I");

  //demux status
  tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b");
  tree->Branch("lowphrat",&lowphrat,"lowphrat/F");
  // David's FD data quality cut
  tree->Branch("dp_fd_quality",&dp_fd_quality,"dp_fd_quality/I");

  //anti numu selection
  float anumupass;
  tree->Branch("anumupass",&anumupass,"anumupass/F");

  //rustem cc pid
  float rustem_cc_pid;
  tree->Branch("rustem_cc_pid",&rustem_cc_pid,"rustem_cc_pid/F");

  //andy cc pid
  float andy_cc_pid;
  tree->Branch("andy_cc_pid",&andy_cc_pid,"andy_cc_pid/F");
  //andy's variables
  float abid_costheta;
  float abid_eventenergy;      
  float abid_trackcharge;   
  float abid_trackenergy;      
  float abid_tracklikeplanes;  
  float abid_trackphfrac;      
  float abid_trackphmean;      
  float abid_trackqpsigmaqp;   
  float abid_y;
  tree->Branch("abid_costheta", &abid_costheta, "abid_costheta/F");
  tree->Branch("abid_eventenergy", &abid_eventenergy, "abid_eventenergy/F");
  tree->Branch("abid_trackcharge", &abid_trackcharge, "abid_trackcharge/I");
  tree->Branch("abid_trackenergy", &abid_trackenergy, "abid_trackenergy/i");
  tree->Branch("abid_tracklikeplanes", &abid_tracklikeplanes, "abid_tracklikeplanes/I");
  tree->Branch("abid_trackphfrac", &abid_trackphfrac, "abid_trackphfrac/F");
  tree->Branch("abid_trackphmean", &abid_trackphmean, "abid_trackphmean/F");
  tree->Branch("abid_trackqpsigmaqp", &abid_trackqpsigmaqp, "abid_trackqpsigmaqp/F");
  tree->Branch("abid_y", &abid_y, "abid_y/F");
                

  float ro_knn_01;
  float ro_knn_10;
  float ro_knn_20;
  float ro_knn_40;
  tree->Branch("ro_knn_01",&ro_knn_01,"ro_knn01/F");
  tree->Branch("ro_knn_10",&ro_knn_10,"ro_knn10/F");
  tree->Branch("ro_knn_20",&ro_knn_20,"ro_knn20/F");
  tree->Branch("ro_knn_40",&ro_knn_40,"ro_knn40/F");

  // beam monitoring variables
  Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill,dt_beam_mon;
  //goodbeam==1 if spill meets criteria defined below 
  //(see the line where goodbeam gets set)
  //beamcnf is a unsigned int as defined in BeamMonSpill::StatusBits
  //  UInt_t beamcnf;
  //now beamcnf is defined by BeamType in conventions
  BeamType::BeamType_t beamcnf;
  //leaving nuTarZ in until we've phased out the summary ntuples
  Double_t nuTarZ;
  //goodbeam==0 if spill doesnt meet these criteria
  Int_t goodbeam;

  //these next variables weren't available in the old beam monitoring ntuples
  UInt_t horn_on, target_in, n_batches;
  Float_t tor101, tr101d, tortgt, trtgtd;
  Float_t hadmon, mumon1, mumon2, mumon3;
  Float_t hposbyspill[6];
  Float_t vposbyspill[6];
  Float_t bibyspill[6];
  int batchnumber;
  tree->Branch("hbw",&hbw,"hbw/D");
  tree->Branch("vbw",&vbw,"vbw/D");
  tree->Branch("hpos1",&hpos1,"hpos1/D");
  tree->Branch("vpos1",&vpos1,"vpos1/D");
  tree->Branch("hpos2",&hpos2,"hpos2/D");
  tree->Branch("vpos2",&vpos2,"vpos2/D");
  tree->Branch("hposbyspill",hposbyspill,"hposbyspill[6]/F");
  tree->Branch("vposbyspill",vposbyspill,"vposbyspill[6]/F");
  tree->Branch("bibyspill",bibyspill,"bibyspill[6]/F");
  tree->Branch("hornI",&hornI,"hornI/D");
  tree->Branch("closestspill",&closestspill,"closestspill/D");
  tree->Branch("dt_beam_mon",&dt_beam_mon,"dt_beam_mon/D");
  tree->Branch("beamcnf",&beamcnf,"beamcnf/i");
  tree->Branch("nuTarZ",&nuTarZ,"nuTarZ/D");
  tree->Branch("goodbeam",&goodbeam,"goodbeam/I");
  tree->Branch("horn_on",&horn_on,"horn_on/i");
  tree->Branch("target_in",&target_in,"target_in/i");
  tree->Branch("n_batches",&n_batches,"n_batches/i");
  tree->Branch("tor101",&tor101,"tor101/F");
  tree->Branch("tr101d",&tr101d,"tr101d/F");
  tree->Branch("tortgt",&tortgt,"tortgt/F");
  tree->Branch("trtgtd",&trtgtd,"trtgtd/F");
  tree->Branch("hadmon",&hadmon,"hadmon/F");
  tree->Branch("mumon1",&mumon1,"mumon1/F");
  tree->Branch("mumon2",&mumon2,"mumon2/F");
  tree->Branch("mumon3",&mumon3,"mumon3/F");
  tree->Branch("batchnumber",&batchnumber,"batchnumber/I");

  //coil status
  Short_t coil_is_ok;
  Short_t coil_is_reverse;
  Short_t coil_stat;
  Float_t coil_current;

  tree->Branch("coil_is_ok",&coil_is_ok,"coil_is_ok/S");
  tree->Branch("coil_is_reverse",&coil_is_reverse,"coil_is_reverse/S");
  tree->Branch("coil_stat",&coil_stat,"coil_stat/S");
  tree->Branch("coil_current",&coil_current,"coil_current/F");

  //gnumiflux info
  Int_t      fluxrun;     // gnumi flux run number
  Int_t      fluxevtno;   // gnumi flux event number
  Float_t    nudxdz;       // dx/dz slope, neutrino rndm decay
  Float_t    nudydz;       // dy/dz slope, neutrino rndm decay
  Float_t    nupz;         // Pz of neutrino (GeV)  rndm decay
  Float_t    nuenergy;     // E(neutrino)    (GeV)  rndm decay
  Float_t    nudxdznear;   // dx/dz slope, neutrino NearDet center
  Float_t    nudydznear;   // dy/dz slope, neutrino NearDet center
  Float_t    nuenergynear;    // E(neutrino)    (GeV)  NearDet center
  Float_t    nwtnear;     // Weight of nu          NearDet center
  Float_t    nudxdzfar;    // dx/dz slope, neutrino FarDet center
  Float_t    nudydzfar;    // dy/dz slope, neutrino FarDet center
  Float_t    nuenergyfar;    // E(neutrino)    (GeV)  FarDet center
  Float_t    nwtfar;      // Weight of nu          FarDet center
  Int_t      norig;       // (ignore)
  Int_t      ndecay;      // Tag of decay mode
  Int_t      ntype;       // Neutrino type (translated to PDG)
  Float_t    vx;          // x vertex of hadron (cm)
  Float_t    vy;          // y vertex of hadron (cm)
  Float_t    vz;          // z vertex of hadron (cm)
  Float_t    pdpx;        // nu parent px at decay point
  Float_t    pdpy;        // nu parent py at decay point
  Float_t    pdpz;        // nu parent pz at decay point
  Float_t    ppdxdz;      // nu parent slope at production
  Float_t    ppdydz;      // nu parent slope at production
  Float_t    pppz;        // nu parent pz at production
  Float_t    ppenergy;    // nu parent energy at production
  Int_t      ppmedium;    // GEANT medium of nu parent at parent production
  Int_t      ptype;       // nu parent type (translated to PDG)
  Float_t    ppvx;        // nu parent production vtx x
  Float_t    ppvy;        // nu parent production vtx y
  Float_t    ppvz;        // nu parent production vtx z
  Float_t    muparpx;     // if parent=mu, hadron parent px
  Float_t    muparpy;     // if parent=mu, hadron parent py
  Float_t    muparpz;     // if parent=mu, hadron parent pz
  Float_t    mupare;      // if parent=mu, hadron parent energy
  Float_t    necm;        // E(nu) in parent cm
  Float_t    nimpwt;      // importance weight
  Float_t    xpoint;      // (unused)
  Float_t    ypoint;      // (unused)
  Float_t    zpoint;      // (unused)
  Float_t    tvx;         // target exit point (x) of parent
  Float_t    tvy;         // target exit point (y) of parent
  Float_t    tvz;         // target exit point (z) of parent
  Float_t    tpx;         // parent px at target exit
  Float_t    tpy;         // parent py at target exit
  Float_t    tpz;         // parent pz at target exit
  Int_t      tptype;      // parent particle type (translated to PDG)
  Int_t      tgen;        // parent generation in cascade

  tree->Branch("fluxrun",&fluxrun,"fluxrun/I");
  tree->Branch("fluxevtno",&fluxevtno,"fluxevtno/I");
  tree->Branch("nudxdz",&nudxdz,"nudxdz/F");
  tree->Branch("nudydz",&nudydz,"nudydz/F");
  tree->Branch("nupz",&nupz,"nupz/F");
  tree->Branch("nuenergy",&nuenergy,"nuenergy/F");
  tree->Branch("nudxdznear",&nudxdznear,"nudxdznear/F");
  tree->Branch("nudydznear",&nudydznear,"nudydznear/F");
  tree->Branch("nuenergynear",&nuenergynear,"nuenergynear/F");
  tree->Branch("nwtnear",&nwtnear,"nwtnear/F");
  tree->Branch("nudxdzfar",&nudxdzfar,"nudxdzfar/F");
  tree->Branch("nudydzfar",&nudydzfar,"nudydzfar/F");
  tree->Branch("nuenergyfar",&nuenergyfar,"nuenergyfar/F");
  tree->Branch("nwtfar",&nwtfar,"nwtfar/F");
  tree->Branch("norig",&norig,"norig/I");
  tree->Branch("ndecay",&ndecay,"ndecay/I");
  tree->Branch("ntype",&ntype,"ntype/I");
  tree->Branch("vx",&vx,"vx/F");
  tree->Branch("vy",&vy,"vy/F");
  tree->Branch("vz",&vz,"vz/F");
  tree->Branch("pdpx",&pdpx,"pdpx/F");
  tree->Branch("pdpy",&pdpy,"pdpy/F");
  tree->Branch("pdpz",&pdpz,"pdpz/F");
  tree->Branch("ppdxdz",&ppdxdz,"ppdxdz/F");
  tree->Branch("ppdydz",&ppdydz,"ppdydz/F");
  tree->Branch("pppz",&pppz,"pppz/F");
  tree->Branch("ppenergy",&ppenergy,"ppenergy/F");
  tree->Branch("ppmedium",&ppmedium,"ppmedium/I");
  tree->Branch("ptype",&ptype,"ptype/I");
  tree->Branch("ppvx",&ppvx,"ppvx/F");
  tree->Branch("ppvy",&ppvy,"ppvy/F");
  tree->Branch("ppvz",&ppvz,"ppvz/F");
  tree->Branch("muparpx",&muparpx,"muparpx/F");
  tree->Branch("muparpy",&muparpy,"muparpy/F");
  tree->Branch("muparpz",&muparpz,"muparpz/F");
  tree->Branch("mupare",&mupare,"mupare/F");
  tree->Branch("necm",&necm,"necm/F");
  tree->Branch("nimpwt",&nimpwt,"nimpwt/F");
  tree->Branch("xpoint",&xpoint,"xpoint/F");
  tree->Branch("ypoint",&ypoint,"ypoint/F");
  tree->Branch("zpoint",&zpoint,"zpoint/F");
  tree->Branch("tvx",&tvx,"tvx/F");
  tree->Branch("tvy",&tvy,"tvy/F");
  tree->Branch("tvz",&tvz,"tvz/F");
  tree->Branch("tpx",&tpx,"tpx/F");
  tree->Branch("tpy",&tpy,"tpy/F");
  tree->Branch("tpz",&tpz,"tpz/F");
  tree->Branch("tptype",&tptype,"tptype/I");
  tree->Branch("tgen",&tgen,"tgen/I");
  
  //pot counting
  float pottor101=0.;
  float pottortgt=0.;
  float pot=0.;
  int NRUNS=0;
  int NSNARLS=0;
  TTree *pottree = new TTree("pottree","pottree");
  pottree->Branch("pot",&pot,"pot/F");
  pottree->Branch("pottor101",&pottor101,"pottor101/F");
  pottree->Branch("pottortgt",&pottortgt,"pottortgt/F");
  pottree->Branch("nruns",&NRUNS,"nruns/I");
  pottree->Branch("nsnarls",&NSNARLS,"nsnarls/I");

  //make a BMSpillAna so we can tell if it's goodbeam
  BMSpillAna bmana;
  bmana.UseDatabaseCuts();

  BeamEnergyCalculator ecalc(&fBECConfig);
  double bec_high=0; double bec_low=0;
  ecalc.GetStandardConfig()->Print();
  ecalc.GetStandardConfig()->Get("beam:high_energy_limit",bec_high);
  ecalc.GetStandardConfig()->Get("beam:low_energy_limit",bec_low);

  // mark's QE id
  //  MadQEID qeid;
  //  qeid.ReadPDFs("/afs/fnal.gov/files/data/minos/d04/kordosky/madqe_files/QEpdfs.root");
  //  bool beamdb=false;
  //  bool checkeddb=false;
  int lastrun=-1;
  for(int ii=0;ii<Nentries;ii++){   
    if(ii%1000==0) std::cout << float(ii)*100./float(Nentries) 
                            << "% done" << std::endl;

    if(!this->GetEntry(ii)) continue;
    NSNARLS++;
    // fill some histograms for mc acceptance
    if(file_is_mc) FillMCHists(&save);

    //give the snarl to the PhysicsNtuple interface
    bool anumuready = phnti->FillSnarl(strecord);
    if(anumuready){
      //      cout<<"ITHINK THE STUPID INTERFACE THING WORKED WITH THE SNARL"<<endl;
    }

    VldContext vc = ntpHeader->GetVldContext();

    nevent = eventSummary->nevent;
    run = ntpHeader->GetRun();
    if(run!=lastrun){
      NRUNS++;
      lastrun = run;
    }
    //get detector type:
    const Detector::Detector_t det 
      = ntpHeader->GetVldContext().GetDetector();


    snarl = ntpHeader->GetSnarl();
    subrun = ntpHeader->GetSubRun();
    vld_sec = ntpHeader->GetVldContext().GetTimeStamp().GetSeconds();
    trgsrc = ntpHeader->GetTrigSrc();
    spilltype = ntpHeader->GetRemoteSpillType();
    trgtime = eventSummary->trigtime;
    litime = eventSummary->litime;
    is_li_sieve = LISieve::IsLI(*strecord);
    coil_stat = detStatus->coilstatus;
    coil_is_ok = 1;
    coil_is_reverse = 0;
    coil_current = 0;
    if(!file_is_mc) {
      coil_is_ok = CoilTools::IsOK(vc);
      coil_is_reverse = CoilTools::IsReverse(vc);
      coil_current = CoilTools::CoilCurrent(vc).first;
    }
    dmx_stat=0;

    //figure out demux status
    dmx_stat+=(dmxStatus->nonphysicalfail);
    dmx_stat+=((dmxStatus->validplanesfail<<1));
    dmx_stat+=((dmxStatus->vertexplanefail<<2));

    //figure out how much of snarl ph was deposited as low ph
    lowphrat=0.;
    lowphrat = ComputeLowPHRatio();

    //is this record in a time period defined as "bad" by Dave Petyt?
    // 1=good, 0=bad    
    if(det==Detector::kFar){
      //      dp_fd_quality=david_passfail(vld_sec);
      //      dp_fd_quality=passfail(vld_sec);
      dp_fd_quality=DataUtil::IsGoodFDData(strecord);
    }
    else if(det==Detector::kNear){
      dp_fd_quality=1;  
    }

    // look through "events" and characterize how close together
    // their vertices are
    // idea is to flag runt events caused by late activity near vertex
    NearbyVertexFinder nby(nevent);
    nby.ProcessEntry(this);
    // run through and tag event with largest signal
    int lgst_evt_idx=-1;
    float big_ph=0.0;
    for(int i=0;i<eventSummary->nevent;i++){ 
      if(!LoadEvent(i)) continue; //no event found
      if(ntpEvent->ph.sigcor > big_ph){ 
        big_ph=ntpEvent->ph.sigcor; 
        lgst_evt_idx=i;
      }
    }

    detector = -1; 
    if(det==Detector::kFar){
      detector = 2;
    }
    else if(det==Detector::kNear){
      detector = 1;     
    }
    // deal with beam type as best we can
    BeamType::BeamType_t current_beam=BeamType::kUnknown;
    if(file_is_mc){
      // for MC set beam to fMCbeam
      current_beam=fMCBeam;
      // for data we will try to deduce from beam monitoring
    }
    
    // first thing, if data do beam monitoring
    hbw=vbw=hpos1=vpos1=hpos2=vpos2=hornI=closestspill=dt_beam_mon=-999.0;
    for(int scnt=0;scnt<6;scnt++){
      hposbyspill[scnt]=0.;
      vposbyspill[scnt]=0.;
      bibyspill[scnt]=0.;
    }
    beamcnf=BeamType::kUnknown;
    nuTarZ=0.0;
    goodbeam=0;
    horn_on=target_in=n_batches=0;
    tor101=tr101d=tortgt=trtgtd=0.0;
    hadmon=mumon1=mumon2=mumon3=0.0;  

    VldTimeStamp vts = vc.GetTimeStamp();

    if(!file_is_mc){
      //      if(beamdb){
      BDSpillAccessor &sa = BDSpillAccessor::Get();
      const BeamMonSpill *bs = sa.LoadSpill(vts);
      hbw=bs->fProfWidX*1000.;//make it in mm
      vbw=bs->fProfWidY*1000.;//make it in mm
      hpos1=bs->fTargProfX*1000.;//make it in mm
      vpos1=bs->fTargProfY*1000.;//make it in mm
      n_batches=0;
      float btint=0.;
      hpos2=0.;
      vpos2=0.;
      for(int bt=0;bt<6;bt++){
        //cout<<"batch "<<bt<<" x "<<bs->fTargBpmX[bt]<<" y "<<bs->fTargBpmY[bt]
        //              <<" int "<<bs->fBpmInt[bt]<<endl;
        hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
        vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);

        hposbyspill[bt]=bs->fTargBpmX[bt];
        vposbyspill[bt]=bs->fTargBpmY[bt];
        bibyspill[bt]=bs->fBpmInt[bt];
        if(bs->fBpmInt[bt]>0.001){
          n_batches++;
        }
        btint+=bs->fBpmInt[bt];
      }
      if(btint!=0){
        hpos2=hpos2*1000./btint;//make it in mm
        vpos2=vpos2*1000./btint;//make it in mm
      }
      else{
        hpos2=-999.;
        vpos2=-999.;
      }
      
      hornI=bs->fHornCur;       
      SpillTimeFinder &sfind= SpillTimeFinder::Instance();

      closestspill = 
        sfind.GetTimeToNearestSpill(vc);
      
      //      beamcnf =bs->GetStatusBits().beam_type;
      beamcnf =bs->BeamType();
      nuTarZ=0.;
      horn_on = bs->GetStatusBits().horn_on;
      target_in = bs->GetStatusBits().target_in;
      
      
      tor101=bs->fTor101;
      tr101d=bs->fTr101d;
      tortgt=bs->fTortgt;
      trtgtd=bs->fTrtgtd;
      hadmon=bs->fHadInt;
      mumon1=bs->fMuInt1;
      mumon2=bs->fMuInt2;
      mumon3=bs->fMuInt3;
      
      goodbeam=0;
      //use Mark D. BMSpillAna to set good beam
      bmana.SetSpill(*bs);
      //      bmana.SetTimeDiff(closestspill);
      // Mark D. suggests that we do the following rather
      // than the line above
      VldTimeStamp delta_vts=vts-bs->SpillTime();
      dt_beam_mon=delta_vts.GetSeconds();
      bmana.SetTimeDiff(delta_vts.GetSeconds());
      
      if(bmana.SelectSpill()){
        goodbeam=1;
      }
      else{
        goodbeam=0;
      }

      if(fDataUseMCBeam==true){
        //allow us to force the data to use the MC beam
        current_beam=fMCBeam;
      }
      else{
        // attempt to deduce beam type from beam monitoring
        //      current_beam = BeamType::FromBeamMon(beamcnf);
        current_beam = beamcnf;
      }
      //      cout<<"Current beam "<<current_beam<<" "<<BeamType::AsString(current_beam)<<endl;

      //don't count fake spills
      if(goodbeam==1&&spilltype!=3){
        //don't count spills where ND coil not working
        if(det==Detector::kFar||coil_is_ok==1){
          pot+=tortgt;
          pottor101+=tor101;
          pottortgt+=tortgt;
        }
      }
    }
    // Write a special record for each spill!
    event=-1; ntrack=-1; nshower=-1;
    reco_emu=0.0; reco_enu=0.0; reco_eshw=0.0; reco_eshw_uncorr=0.0;
    true_enu=0.0; true_emu=0.0; true_eshw=0.0;
    pass=0; is_fid=0; is_fid_2007=0;
    anumupass=0;
    rustem_cc_pid=0;


    ro_knn_01=0.;
    ro_knn_10=0.;
    ro_knn_20=0.;
    ro_knn_40=0.;

    andy_cc_pid=-1.;
    abid_costheta         =0.;
    abid_eventenergy      =0.;
    abid_trackcharge      =0;
    abid_trackenergy      =0;
    abid_tracklikeplanes  =0;
    abid_trackphfrac      =0.;
    abid_trackphmean      =0.;
    abid_trackqpsigmaqp   =0.;
    abid_y                =0.;

    tree->Fill();

    if(!openedabpidfile){
      string base="";
      base=getenv("SRT_PRIVATE_CONTEXT");
      if(base!=""&&base!="."){
        // check if directory exists in SRT_PRIVATE_CONTEXT      
        std::string path = base + "/Mad/data";
        void *dir_ptr = gSystem -> OpenDirectory(path.c_str());
        if(!dir_ptr){
          // if it doesn't exist then use SRT_PUBLIC_CONTEXT
          base=getenv("SRT_PUBLIC_CONTEXT"); 
        }
      }
      else{
        base=getenv("SRT_PUBLIC_CONTEXT");
      }
      if(base=="") {
        cout<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
        assert(false);
      }
      base+="/Mad/data";
      
      string fname=base;
      string sdet="";
      if(det==Detector::kNear){
        sdet="near";
      }
      else if(det==Detector::kFar){
        sdet="far";
      }
      else{
        cout<<"Could not set ab pid file, defaulting to far"<<endl;
        sdet="far";
      }      
      string rmc="";
      if(ReleaseType::IsCedar(reltype)&&mcver=="daikon"){      
        rmc="cedar_daikon";
      }
      else{
        cout<<"Dont know reco/mc versions "
            <<"defaulting to cedar_daikon"<<endl;
        rmc="cedar_daikon";
      }
      string sbeam="";
      switch (current_beam){
      case BeamType::kL010z000i:
        sbeam="le0";
        break;
      case BeamType::kL010z170i:
        sbeam="le170";
        break;
      case BeamType::kL010z185i:
        sbeam="le";
        break;
      case BeamType::kL010z200i:
        sbeam="le200";
        break;
      case BeamType::kL100z200i:
        sbeam="pme";
        break;
      case BeamType::kL150z200i:
        sbeam="pmhe";
        break;
      case BeamType::kL250z200i:
        sbeam="phe";
        break;
      case BeamType::kLE:
        sbeam="le";
        break;
      default:
        cout<<"Don't know beam type "
            <<" defaulting to LE"<<endl;
        sbeam="le";
        break;
      }
      fname+="/ab_pdf_"+sdet+"_"+sbeam+"_"+rmc+".root";
      cout<<"Reading AB PDFS from "<<fname<<endl;
      abid.ReadPDFs(fname.c_str());
      openedabpidfile=true;
    }


    if(nevent==0) continue;    
    //fill tree once for each reconstructed event
    for(int i=0;i<eventSummary->nevent;i++){ 
      if(!LoadEvent(i)) continue; //no event found

      //set event index:
      event = i;
      ntrack = ntpEvent->ntrack;
      nshower = ntpEvent->nshower;
      nstrip = ntpEvent->nstrip;

      //anti nu selection
      anumupass=0;
      rustem_cc_pid=0;
      ro_knn_01=0.;
      ro_knn_10=0.;
      ro_knn_20=0.;
      ro_knn_40=0.;

      andy_cc_pid=-1.;
      abid_costheta         =0.;
      abid_eventenergy      =0.;
      abid_trackcharge      =0;
      abid_trackenergy      =0;
      abid_tracklikeplanes  =0;
      abid_trackphfrac      =0.;
      abid_trackphmean      =0.;
      abid_trackqpsigmaqp   =0.;
      abid_y                =0.;
      
      if(anumuready){
        //      cout<<"Filling event rustem stuff"<<endl;

        anumupass     = phnti->GetVar("numubar", ntpEvent);
        rustem_cc_pid = phnti->GetVar("knn_pid", ntpEvent);
        ro_knn_01     = phnti->GetVar("knn_01",  ntpEvent);
        ro_knn_10     = phnti->GetVar("knn_10",  ntpEvent);
        ro_knn_20     = phnti->GetVar("knn_20",  ntpEvent);
        ro_knn_40     = phnti->GetVar("knn_40",  ntpEvent);

      }


      
      //zero all tree values
      true_enu = 0; true_emu = 0; true_eshw = 0; true_x = 0; true_y = 0; 
      true_q2 = 0; true_w2 = 0; true_dircosneu = 0; true_dircosz = 0;
      true_vtxx = 0; true_vtxy = 0; true_vtxz = 0; true_pitt_fid=0;
      true_trk_cmplt=0.;
      resnum=0;
      npp=npm=npz=np=nn=0;
      epp=epm=epz=ep=en=0.0;

      fluxrun=0; fluxevtno=0; nudxdz=0.; nudydz=0.; nupz=0.; nuenergy=0.;
      nudxdznear=0.; nudydznear=0.; nuenergynear=0.; nwtnear=0.; nudxdzfar=0.;
      nudydzfar=0.; nuenergyfar=0.; nwtfar=0.; norig=0; ndecay=0; ntype=0;
      vx=0.; vy=0.; vz=0.; pdpx=0.; pdpy=0.; pdpz=0.; ppdxdz=0.; ppdydz=0.; pppz=0.;
      ppenergy=0.; ppmedium=0; ptype=0; ppvx=0.; ppvy=0.; ppvz=0.;
      muparpx=0.; muparpy=0.; muparpz=0.; mupare=0.; necm=0.; nimpwt=0.;
      xpoint=0.; ypoint=0.; zpoint=0.; tvx=0.; tvy=0.; tvz=0.; tpx=0.;
      tpy=0.; tpz=0.; tptype=0; tgen=0;
      
      flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0;
      rc_boundary=0;
      mcevent = -1;
      is_fid = 0;is_fid_2007 = 0; is_cev = 0; sr_contained = 0; is_mc = 0; pass_basic = 0; 
      is_pitt_fid=0; pitt_evt_class=0; pitt_trk_qual=0; nd_radial_fid=0;
      is_trk_fid=0;
      duvvtx=0;
      pid0 = 0; pid1 = 0; pid2 = 0; pass = 0;
      reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_wtd_eshw=0; 
      reco_eshw_uncorr=0; reco_wtd_eshw_uncorr=0;
      reco_vtx_eshw=0; reco_vtx_wtd_eshw=0;
      reco_qe_enu = 0;  reco_mushw=0;
      reco_qe_q2 = 0; reco_dircosneu = 0;
      reco_dircosx = 0;; reco_dircosy = 0;; reco_dircosz = 0;
      reco_vtxx = reco_vtxy = reco_vtxz = reco_vtxr=0; 
      reco_trk_vtxx = reco_trk_vtxy = reco_trk_vtxz = reco_trk_vtxr=0; 
      reco_shw_vtxx = reco_shw_vtxy = reco_shw_vtxz = reco_shw_vtxr=0; 
      reco_trk_endx = reco_trk_endy = reco_trk_endz = reco_trk_endr=0; 
      evt_nstrips = slc_nstrips = trk_nstrips=0;
      evt_nplanes = slc_nplanes = trk_nplanes=0;
      shw_nplanes_cut=shw_nplanes_raw=0;
      shw_ph_cut=shw_ph_raw=0.0;

      evtlength = trklength = shwlength=trknstp=0; 
      trkphfrac = trkphplane = 0;
      ntrklike= trkend=trkendu=trkendv=shwbeg=shwend=0;
      trkmomrange = 0; trkqp = 0; trkeqp = 0;
      trkchi2=0.0; trkndf=0;
      pass_track=0; emu_meth=0;
      trk_fit_pass =0;
      // kinematic quantities always positive
      // convention: -1 indicates unfilled variable (NC events)
      reco_y = reco_x = reco_Q2 = reco_W2 = -1;
      
      // vertex back tracking

      nvsi=nvti=bvsi=bvti=-1;
      nvsevt=bvsevt=nvtevt=bvtevt=-1;
      nvsd=nvtd=bvsd=bvtd=999999;
      nvst=bvst=nvtt=bvtt=1e6;
      nvsp=nvtp=bvsp=bvtp=-1.0;

      evtph=0; evtphgev=0;
      evtsigfull=evtsigpart=trksigfull=trksigpart=0;            

      evttimemin=evttimemax=evttimeln=0.;
      evttimeminphc=evttimemaxphc=evttimelnphc=0.;

      earlywph=earlywphpw=earlywphsphere=0.;
      ph1000=ph500=ph100=0.;
      lateph1000=lateph500=lateph100=lateph50=0.;


      nduprs=0;
      dupphrs=0.;
      
      nu_px=nu_py=nu_pz=tar_px=tar_py=tar_pz=tar_e=0.0;
      inu=0;
      inunoosc=0;
      ntotss=0;
      for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;}

      med_qe_pid=-1001.0;
      med_qe_pass=0;
      niki_cc_pid=-999;
      dave_cc_pid=-999;
            
      //check for a corresponding mc event
      // first, looks to match reconstructed track
      // then, looks to match reconstructed shower
      if(LoadTruthForRecoTH(i,mcevent)) {
        
        // should flag with a "bad truth word"
        // in case the function above returns false

        //true info for tree:
        is_mc             = 1;
        nucleus           = Nucleus(mcevent);
        flavour           = Flavour(mcevent);
        initial_state     = Initial_state(mcevent);
        process           = IResonance(mcevent); 
        cc_nc             = IAction(mcevent);
        true_enu          = TrueNuEnergy(mcevent); 
        true_emu          = TrueMuEnergy(mcevent); 
        true_eshw         = TrueShwEnergy(mcevent); 
        true_x            = X(mcevent);
        true_y            = Y(mcevent);
        true_q2           = Q2(mcevent);
        true_w2           = W2(mcevent);
        true_dircosz      = TrueMuDCosZ(mcevent);
        true_dircosneu    = TrueMuDCosNeu(mcevent);
        

        //flux info
        if(LoadTruth(mcevent)){
          fluxrun=ntpTruth->flux.fluxrun;
          fluxevtno=ntpTruth->flux.fluxevtno;
          nudxdz=ntpTruth->flux.ndxdz;
          nudydz=ntpTruth->flux.ndydz;
          nupz=ntpTruth->flux.npz;
          nuenergy=ntpTruth->flux.nenergy;
          nudxdznear=ntpTruth->flux.ndxdznear;
          nudydznear=ntpTruth->flux.ndydznear;
          nuenergynear=ntpTruth->flux.nenergynear;
          nwtnear=ntpTruth->flux.nwtnear;
          nudxdzfar=ntpTruth->flux.ndxdzfar;
          nudydzfar=ntpTruth->flux.ndydzfar;
          nuenergyfar=ntpTruth->flux.nenergyfar;
          nwtfar=ntpTruth->flux.nwtfar;
          norig=ntpTruth->flux.norig;
          ndecay=ntpTruth->flux.ndecay;
          ntype=ntpTruth->flux.ntype;
          vx=ntpTruth->flux.vx;
          vy=ntpTruth->flux.vy;
          vz=ntpTruth->flux.vz;
          pdpx=ntpTruth->flux.pdpx;
          pdpy=ntpTruth->flux.pdpy;
          pdpz=ntpTruth->flux.pdpz;
          ppdxdz=ntpTruth->flux.ppdxdz;
          ppdydz=ntpTruth->flux.ppdydz;
          pppz=ntpTruth->flux.pppz;
          ppenergy=ntpTruth->flux.ppenergy;
          ppmedium=ntpTruth->flux.ppmedium;
          ptype=ntpTruth->flux.ptype;
          ppvx=ntpTruth->flux.ppvx;
          ppvy=ntpTruth->flux.ppvy;
          ppvz=ntpTruth->flux.ppvz;
          muparpx=ntpTruth->flux.muparpx;
          muparpy=ntpTruth->flux.muparpy;
          muparpz=ntpTruth->flux.muparpz;
          mupare=ntpTruth->flux.mupare;
          necm=ntpTruth->flux.necm;
          nimpwt=ntpTruth->flux.nimpwt;
          xpoint=ntpTruth->flux.xpoint;
          ypoint=ntpTruth->flux.ypoint;
          zpoint=ntpTruth->flux.zpoint;
          tvx=ntpTruth->flux.tvx;
          tvy=ntpTruth->flux.tvy;
          tvz=ntpTruth->flux.tvz;
          tpx=ntpTruth->flux.tpx;
          tpy=ntpTruth->flux.tpy;
          tpz=ntpTruth->flux.tpz;
          tptype=ntpTruth->flux.tptype;
          tgen=ntpTruth->flux.tgen;       
        }

        if(ntpTHTrack){
          if(ntpTHTrack->neumc==mcevent){
            true_trk_cmplt=ntpTHTrack->completeslc;
          }
        }
        
        Float_t* true_p = TrueNuMom(mcevent);
        if(true_p){
          nu_px=true_p[0]; nu_py=true_p[1]; nu_pz=true_p[2];
          delete[] true_p; true_p=0;
        }
        true_p = Target4Vector(mcevent);
        if(true_p){
          tar_px=true_p[0]; tar_py=true_p[1]; tar_pz=true_p[2]; 
          tar_e=true_p[3];
          delete[] true_p; true_p=0;
        }

        Float_t *vtx = TrueVtx(mcevent);
        true_vtxx         = vtx[0];
        true_vtxy         = vtx[1];
        true_vtxz         = vtx[2];
        float true_vtxu = (true_vtxx + true_vtxy)/sqrt(2.0);
        float true_vtxv = (true_vtxy - true_vtxx)/sqrt(2.0);
        true_pitt_fid=PittInFidVol(true_vtxx, true_vtxy, 
                                   true_vtxz,
                                   true_vtxu, true_vtxv);

        resnum = HadronicFinalState(mcevent);
        npp = NumFSPion(mcevent,+1);
        npm = NumFSPion(mcevent,-1);
        npz = NumFSPiZero(mcevent);
        np = NumFSProton(mcevent);
        nn = NumFSNeutron(mcevent);
        
        epp = TotFSPionEnergy(mcevent,+1);
        epm = TotFSPionEnergy(mcevent,-1);
        epz = TotFSPiZeroEnergy(mcevent);
        ep = TotFSProtonEnergy(mcevent);
        en = TotFSNeutronEnergy(mcevent);
        
        // beam energy weights //
        inu=0;
        inu = INu(mcevent);
        inunoosc = INuNoOsc(mcevent);
        bec_inve=bec_le=bec_050=bec_100=bec_200=bec_250=1.0;
        /*
        bec_inve=ecalc.GetWeight(BeamType::kInverseE,det,inu,
                                 true_enu, bec_high,bec_low);
        bec_le=ecalc.GetWeight(BeamType::kLE,det,inu,
                               true_enu, bec_high,bec_low);
        bec_050=ecalc.GetWeight(BeamType::k050,det,inu,
                                 true_enu, bec_high,bec_low);
        bec_100=ecalc.GetWeight(BeamType::k100,det,inu,
                                 true_enu, bec_high,bec_low);
        bec_200=ecalc.GetWeight(BeamType::k200,det,inu,
                                 true_enu, bec_high,bec_low);
        bec_250=ecalc.GetWeight(BeamType::k250,det,inu,
                                 true_enu, bec_high,bec_low);
        */
        delete [] vtx;
      }

      // is the event in the fiducial volume
      // MadQuantities::IsFidVtxEvt()
      //  this is used by Dave/Niki analysis
      // provisional, will update to using track vertices
      // if there is a track.      
      is_fid_2007 = IsFidVtxEvt(event);
      is_fid = IsFid_2008(event);
      if(det==Detector::kFar){
        if(IsFidFD_AB(event)){
          is_fid_2007 = 1;
        }
        else{
          is_fid_2007=0;
        }
      }

      // is this event on a readout crate boundary
      rc_boundary=0;      
      if(det==Detector::kFar){
        //      if(FDRCBoundary(eventSummary->planeall.beg,eventSummary->planeall.end)) rc_boundary=1;
        //      if(FDRCBoundary(eventSummary->plane.beg,eventSummary->plane.end)) rc_boundary=1;
        rc_boundary=FDRCBoundary();
      }
      // is this event the largest (sicor) in the spill
      if(event==lgst_evt_idx) largest_event=1;
      else largest_event=0;
      

      if(PassBasicCuts(event)) { // this does nothing.
        //reco info for tree:
        pass_basic        = 1;
        
        
        //index will -1 if no track/shower present in event
        int track_index   = -1;
        // track spanning largest number of planes
        LoadLargestTrackFromEvent(i,track_index);
        if(track_index>-1&&det==Detector::kNear){
          //continue to use PRL cuts for ND fid. volume
          //FD fid vol has been updated and is computed above
          is_fid_2007=InFidVol(ntpTrack->vtx.x,ntpTrack->vtx.y,ntpTrack->vtx.z);
        }
        pass_track = PassTrackCuts(track_index);
        //methods should all be protected against index = -1
        
        //      reco_emu          = RecoMKMuEnergy(emu_meth,track_index,
        //                                         !file_is_mc,det);
        reco_emu = RecoMuEnergyNew(0,track_index,vc,
                                   reltype,EnergyCorrections::kDefault);

        int shower_index  = -1;
        

        //      LoadLargestShowerFromEvent(i,shower_index);
        // Now use 0th shower (set by Jim) with cuts to remove brems
        LoadShower_Jim(i,shower_index,det);
        // shower with the best match to the vertex
        //      reco_eshw  = RecoShwEnergy(shower_index, det,1,!file_is_mc);

        reco_eshw=RecoShwEnergyNew(shower_index,0,vc,
                                    reltype,EnergyCorrections::kDefault);

        reco_wtd_eshw=RecoShwEnergyNew(shower_index,1,vc,
                                       reltype,EnergyCorrections::kDefault);
        
        //      cout<<"reco_eshw output"<<reco_eshw<<endl;
        //      cout<<"VLD context "<<endl;
        //      vc.Print();
        //      cout<<endl;


         if(shower_index>-1) {
           reco_eshw_uncorr = ntpShower->shwph.linCCgev;
           reco_wtd_eshw_uncorr = ntpShower->shwph.wtCCgev;
         }
         reco_vtx_eshw = reco_eshw;
         reco_vtx_wtd_eshw = reco_wtd_eshw;

         // provisional, to be updated later
         reco_enu = reco_emu + reco_eshw;
         // event energy according to offline
         evtphgev = ntpEvent->ph.gev;
         if(reco_enu>0){
           reco_qe_enu       = RecoQENuEnergy(track_index);
           reco_qe_q2        = RecoQEQ2(track_index);
           // note, NtpSRFiducial isn't filled for ND events
           // is_cev determines if the track is a stopper or punch-through
           // used in Dave/Niki analysis
           is_cev            = IsFidAll(track_index);
           // use Pitt group's fiducial cuts for ND events
           if(detector==1){
             // for now, use event vertex
             // but maybe should use track vertex...
             // if vertex finder works well it shouldn't matter...
             is_pitt_fid = PittInFidVol(ntpEvent->vtx.x, ntpEvent->vtx.y, 
                                        ntpEvent->vtx.z,
                                        ntpEvent->vtx.u, ntpEvent->vtx.v);
             nd_radial_fid = NDRadialFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,
                                            ntpEvent->vtx.z);
           }
           else{
             is_pitt_fid = InFidVol();
           }
           // have already checked that ntpEvent exists
           reco_vtxx         = ntpEvent->vtx.x;
           reco_vtxy         = ntpEvent->vtx.y;
           reco_vtxz         = ntpEvent->vtx.z;

           // check if distance from vertex to shower is too large
           // if so, set reco_eshw =0 
           // idea is to handle case of no vertex shower but runt downstream
           if(shower_index!=-1){
             //     bool shw_ok=true;
             const float max_shw_dist=14*0.06; // about 2 interaction lengths
             float dist_z = fabs(reco_vtxz - ntpShower->vtx.z);
             float dist_r = sqrt( pow(ntpEvent->vtx.u-ntpShower->vtx.u,2) +
                                  pow(ntpEvent->vtx.v-ntpShower->vtx.v,2) );
             if( (dist_z + dist_r) > max_shw_dist ){
               reco_vtx_eshw=0;
               reco_vtx_wtd_eshw=0;
             }

           }
           // radial position of vertex
           if(detector==1){
             reco_vtxr = pow(reco_vtxx-kXcenter,2) 
               + pow(reco_vtxy-kYcenter,2);
             reco_vtxr=sqrt(reco_vtxr);
           }
           else reco_vtxr = sqrt (pow(reco_vtxx,2)+pow(reco_vtxy,2));

           // compute signal in fully and partially instrumented regions
           SigInOut(track_index,evtsigfull,evtsigpart,trksigfull,trksigpart);
            if(detector==2){
             //if far det, all signal is in fully instrumented region
             evtsigfull=ntpEvent->ph.sigcor;
             evtsigpart=0.;
             if(track_index!=-1){
               trksigfull=ntpTrack->ph.sigcor;
             }
             else{
               trksigfull=0.;
             }
             trksigpart=0.;
           }

           // angle reconstruction
           if(track_index>-1) {
             reco_dircosx = ntpTrack->vtx.dcosx;
             reco_dircosy = ntpTrack->vtx.dcosy;
           }
           reco_dircosz      = RecoMuDCosZ(track_index);

           // event vertex
           float vtx_array[3]; 
           vtx_array[0]=ntpEvent->vtx.x; vtx_array[2]=ntpEvent->vtx.y;
           vtx_array[2]=ntpEvent->vtx.z;          
           if(detector==1) 
             reco_dircosneu = RecoMuDCosNeuND(track_index, vtx_array);
           else reco_dircosneu = RecoMuDCosNeuFD(track_index, vtx_array);

           evtlength         = ntpEvent->plane.end-ntpEvent->plane.beg;
           if(track_index!=-1){ //check track is present before using ntpTrack

             trklength         = ntpTrack->plane.end-ntpTrack->plane.beg;
             ntrklike =ntpTrack->plane.ntrklike;
             trkend         = ntpTrack->plane.end;
             trkendu         = ntpTrack->plane.endu;
             trkendv         = ntpTrack->plane.endv;
             trknstp       = ntpTrack->nstrip;
             trkmomrange       = ntpTrack->momentum.range;

             trk_fit_pass       = ntpTrack->fit.pass;
             //     trkmomrange = EnergyCorrections::CorrectMomentumFromRange(trkmomrange,!file_is_mc,det);
             //now use FullyCorrect Versions
             trkmomrange = EnergyCorrections::FullyCorrectMomentumFromRange(trkmomrange,vc,reltype);
             trkqp             = ntpTrack->momentum.qp;
             if(trkqp!=0){
               //             trkqp = 1.0/EnergyCorrections::CorrectSignedMomentumFromCurvature(1.0/trkqp,
               trkqp = 1.0/EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(1.0/trkqp,vc,reltype);
             }

             trkeqp            = ntpTrack->momentum.eqp;
             trkphfrac         = ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
             trkphplane        = ntpTrack->ph.mip/ntpTrack->plane.n;
             trkchi2           = ntpTrack->fit.chi2;
             trkndf           = ntpTrack->fit.ndof;
             reco_trk_vtxx         = ntpTrack->vtx.x;
             reco_trk_vtxy         = ntpTrack->vtx.y;
             reco_trk_vtxz         = ntpTrack->vtx.z;

             reco_trk_endx         = ntpTrack->end.x;
             reco_trk_endy         = ntpTrack->end.y;
             reco_trk_endz         = ntpTrack->end.z;

             sr_contained          = ntpTrack->contained;
             
             // compute CC event class as defined by Pitt group
             pitt_evt_class = PittTrkContained(ntpTrack->end.x,
                                               ntpTrack->end.y,
                                               ntpTrack->end.z,
                                               ntpTrack->end.u,
                                               ntpTrack->end.v);
             duvvtx=abs(ntpTrack->plane.begu - ntpTrack->plane.begv);
             // based on the event class, go back and recompute
             // the muon and neutrino energy
             // modified here, now use is_cev=IsFidAll()
             int do_meth=0; // emu reco method we should use
             if( is_cev==1 ) 
               do_meth=2; // stoppers --> range
             else if( is_cev==0 ) 
               do_meth=1; // punch-through --> curvature

             // with the advent of trk.fit.pass==0 being accepted
             // it is possible that we have to modify do_meth to
             // account for qp==0
             // in this case we have to use the range
             if(trkqp==0) do_meth=2;          

             if(do_meth==1 || do_meth==2){
               //             reco_emu = RecoMKMuEnergy(do_meth,track_index,!file_is_mc,det);
               reco_emu = RecoMuEnergyNew(do_meth,track_index,
                                          vc,reltype,
                                          EnergyCorrections::kDefault);
               // update reco_enu to account for possible change 
               // in reco_emu
               reco_enu = reco_eshw + reco_emu;
             }
             emu_meth=do_meth; // was forgetting about this

             //check to see if track vertex is in fid. vol
             if(detector==1){
               is_trk_fid = PittInFidVol(ntpTrack->vtx.x, ntpTrack->vtx.y, 
                                         ntpTrack->vtx.z,
                                         ntpTrack->vtx.u, ntpTrack->vtx.v);
             }
             else{
               is_trk_fid = InFidVol(ntpTrack->vtx.x,ntpTrack->vtx.y,
                                     ntpTrack->vtx.z);
             }

             // look for showers along muon track
             reco_mushw=MuonShowerEnergy(ntpEvent, ntpTrack);

             // pitt tracking quality cuts
             Int_t temp_ndof = ntpTrack->fit.ndof;
             if(temp_ndof==0) temp_ndof=1;
             if( (ntpTrack->fit.pass>0) 
                 && (ntpTrack->fit.chi2/temp_ndof <20 )
                 && abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
               pitt_trk_qual=1;

             const double M=(0.93827 + 0.93957)/2.0; // <nucleon mass>
             if(reco_emu>0) reco_Q2 = 
                              2*reco_enu*reco_emu*(1.0 - reco_dircosneu);
             reco_W2 = M*M - reco_Q2 + 2*M*reco_eshw;     
             if(reco_enu>0) reco_y = reco_eshw/reco_enu;
             if(reco_eshw>0 && reco_Q2>0) reco_x =  reco_Q2/(2*M*reco_eshw);

             const double muM=0.10566; // muon mass
             const double muP=sqrt(reco_emu*reco_emu -muM*muM);
             reco_qe_enu = (M*reco_emu - muM*muM/2.0)
               /(M - reco_emu + muP*reco_dircosneu);
             reco_qe_q2 = abs(-2.0*reco_qe_enu*(reco_emu-muP*reco_dircosneu)+muM*muM);


             med_qe_pid=qeid.AltCalcQEID(this,event,reco_enu,reco_W2);
             med_qe_pass=qeid.PassMEDQECut(reco_enu,med_qe_pid);

             if(nsid.ChooseWeightFile(det,current_beam)){
               if(!nsid.GetPID(this,event,det,niki_cc_pid)) niki_cc_pid=-999;
             }
             else niki_cc_pid=-999;
             // original code made detector dependent fiducial volume cut here
             // but, it's not necessary since PID has already been
             // calculated...
             // must reevaluate efficiency or make same cut downstream      

             if(det==Detector::kFar && !file_is_mc){
               if(recover==""||recover.find("birch")<recover.size()||
                  recover.find("Birch")<recover.size()||
                  recover.find("BIRCH")<recover.size()||
                  recover.find("R1.18")<recover.size()){
                 dpid.SetPHCorrection(1.018);
               }
               if(recover.find("cedar")<recover.size()||
                  recover.find("Cedar")<recover.size()||
                  recover.find("CEDAR")<recover.size()||
                  recover.find("R1.24")<recover.size()){
                 dpid.SetPHCorrection(1.0);
               }
             }
             if(dpid.ChoosePDFs(det,current_beam,recover,mcver)){
               dave_cc_pid = dpid.CalcPID(this,event,0);
             }


             //     if(abid.PassCuts(ntpEvent,strecord)){
             andy_cc_pid=abid.CalcPID(ntpEvent,strecord);
             abid_costheta         =abid.GetRecoCosTheta(ntpEvent,strecord);
             abid_eventenergy      =abid.GetRecoE(ntpEvent,strecord);
             abid_trackcharge      =abid.GetTrackEMCharge(ntpEvent,strecord);
             abid_trackenergy      =abid.GetTrackPlanes(ntpEvent,strecord);
             abid_tracklikeplanes  =abid.GetTrackLikePlanes(ntpEvent,strecord);
             abid_trackphfrac      =abid.GetTrackPHfrac(ntpEvent,strecord);
             abid_trackphmean      =abid.GetTrackPHmean(ntpEvent,strecord);
             abid_trackqpsigmaqp   =abid.GetTrackQPsigmaQP(ntpEvent,strecord);
             abid_y                =abid.GetRecoY(ntpEvent,strecord);
             //     }

             TObject *recobj=0;
             if(record!=0){
               recobj=record;
             }
             else{
               recobj=strecord;
             }


             LoadTrack(track_index);
             //figure out planes in trk > 1.5 pe
             float temp_ph;
             trk_nplanes=NPlanesInObj(recobj,ntpTrack,1.5,temp_ph);
             //figure out strips in trk > 1.5 pe
             trk_nstrips=NStripsInObj(recobj,ntpTrack,1.5,temp_ph);

           }
           else {
             trklength         = 0;
             trkmomrange       = 0;
             trkqp             = 0;
             trkeqp            = 0;
             trkphfrac         = 0;
             trkphplane        = 0;
             is_trk_fid        = -1;
           }

           //this used to load the slice associated with the track
           //but we really want the slice associated witht the event
           //tv 10-19-05
           LoadSlice(ntpEvent->slc);
           //       cout<<"*********SNARL="<<snarl<<" "<<event
           //           <<"*********************************"<<endl;
           TObject *recobj=0;
           if(record!=0){
             recobj=record;
           }
           else{
             recobj=strecord;
           }

           float temp_ph;
           //figure out planes in slice > 1.5 pe
           slc_nplanes=NPlanesInObj(recobj,ntpSlice,1.5,temp_ph);
           //figure out planes in evt > 1.5 pe
           evt_nplanes=NPlanesInObj(recobj,ntpEvent,1.5,temp_ph);
           //figure out strips in slice > 1.5 pe
           slc_nstrips=NStripsInObj(recobj,ntpSlice,1.5,temp_ph);
           //figure out strips in evt > 1.5 pe
           evt_nstrips=NStripsInObj(recobj,ntpEvent,1.5,temp_ph);


           if(shower_index!=-1 && reco_eshw>0){ // case of no vertex shower
             LoadShower(shower_index);
             shwlength=ntpShower->plane.n;
             shw_nplanes_raw=NPlanesInObj(recobj,ntpShower,0.0,shw_ph_raw);
             shw_nplanes_cut=NPlanesInObj(recobj,ntpShower,1.5,shw_ph_cut);

             shwend=ntpShower->plane.end;
             shwbeg=ntpShower->plane.beg;
             reco_shw_vtxx=ntpShower->vtx.x;
             reco_shw_vtxy=ntpShower->vtx.y;
             reco_shw_vtxz=ntpShower->vtx.z;


             int ncl = ntpShower->ncluster;

             //Float_t  zvtx
             //Float_t  tposvtx
             //Float_t  slope
             //Float_t  avgdev
             /*
             static TH1F hswu("hswu","",400,4,-4);
             static TH1F hswu_wtd("hswu_wtd","",400,4,-4);
             static TH1F hswv("hswv","",400,4,-4);
             static TH1F hswv_wtd("hswv_wtd","",400,4,-4);
             hswu.Reset();
             hswu_wtd.Reset();
             hswv.Reset();
             hswv_wtd.Reset();
             */
             Moments mom_swu;
             Moments mom_swu_wtd;
             Moments mom_swv;
             Moments mom_swv_wtd;

             etu=etv=elu=elv=0.0; 
             swu=swv=wswu=wswv=0.0;

             //     const float ss_cut=0.150;//GeV
             const float ss_cut=0.0;//GeV
             float totsigssu=0.0; // em and hadr subshowers
             float totsigssv=0.0; // em and hadr subshowers
             for(int ii=0; ii<ncl; ii++){
               if(LoadCluster(ntpShower->clu[ii])){
                 if(ntpCluster->ph.gev >ss_cut){
                   nss[ntpCluster->id]+=1;
                   ess[ntpCluster->id]+=ntpCluster->ph.gev;
                   ntotss++;
                   if(ntpCluster->id <2){
                     // slope = tpos/zpos?
                     float sinang = sin(atan(ntpCluster->slope));
                     if(ntpCluster->planeview==PlaneView::kU){
                       totsigssu+=ntpCluster->ph.gev;
                       //hswu.Fill(ntpCluster->tposvtx);
                       mom_swu.AddPoint(ntpCluster->tposvtx);
                       mom_swu_wtd.AddPoint(ntpCluster->tposvtx,
                                            ntpCluster->ph.gev);

                       etu += sinang*ntpCluster->ph.gev;
                       elu += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev;
                     }
                     else if(ntpCluster->planeview==PlaneView::kV){
                       totsigssv+=ntpCluster->ph.gev;
                       //hswv.Fill(ntpCluster->tposvtx);
                       mom_swv.AddPoint(ntpCluster->tposvtx);
                       mom_swv_wtd.AddPoint(ntpCluster->tposvtx,
                                            ntpCluster->ph.gev);

                       etv += sinang*ntpCluster->ph.gev;
                       elv += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev;
                     }
                   }
                 }
               }
             }
             /*
             for(int ii=0; ii<ncl; ii++){
               if(LoadCluster(ntpShower->clu[ii])){
                 if(ntpCluster->ph.gev >ss_cut){
                   if(ntpCluster->id <2){
                     if(ntpCluster->planeview==PlaneView::kU){
                       totsigssu+=ntpCluster->ph.gev;
                       hswu_wtd.Fill(ntpCluster->tposvtx);
                     }
                     else if(ntpCluster->planeview==PlaneView::kV){
                       totsigssv+=ntpCluster->ph.gev;
                       hswv_wtd.Fill(ntpCluster->tposvtx);
                     }
                   }
                 }
               }
             }
             */
             /*
             swv=hswv.GetRMS();
             swu=hswu.GetRMS();
             wswv=hswv_wtd.GetRMS();
             wswu=hswu_wtd.GetRMS();
             */
             swv=mom_swv.GetRMS();
             swu=mom_swu.GetRMS();
             wswv=mom_swv_wtd.GetRMS();
             wswu=mom_swu_wtd.GetRMS();

           }
           /*
           if(ntrack==1 && reco_enu>0.0 && reco_enu<140
              && is_pitt_fid && pitt_trk_qual && (is_trk_fid!=0)) pass=1;
           */
           if(ntrack>0 && trk_fit_pass==1 && is_fid==1) pass=1;


         } // if(reco_enu>0)
       } // if(PassBasicCuts())

       // for both near in time and near in space
       // index of nearest, dist of nearest, ph of nearest
       // near_vtx_space_idx, near_vtx_space_dist, near_vtx_space_ph 
       // nvsi, nvsd, nvsp
       // near_vtx_time_idx, near_vtx_time_dist, near_vtx_time_ph 
       // nvti, nvtd, nvtp
       // for the nearest event, the index of the event nearest it
       // ph of event nearest it, distance of event nearest it
       //  maybe use for event rehashing
       // back_vtx_space_idx, back_vtx_space_dist, back_vtx_space_ph 
       // bvsi, bvsd, bvsp
       // back_vtx_time_idx, back_vtx_time_dist, back_vtx_time_ph 
       // bvti, bvtd, bvtp

       nby.ClosestSpatial(i,nvsi,nvsd,nvst,nvsp,nvsevt);
       nby.ClosestSpatial(nvsi,bvsi,bvsd,bvst,bvsp,bvsevt);
       nby.ClosestTemporal(i,nvti,nvtd,nvtt,nvtp,nvtevt);
       nby.ClosestTemporal(nvti,bvti,bvtd,bvtt,bvtp,bvtevt);

       GetEvtTimeChar(evttimemin,evttimemax,evttimeln);
       GetEvtTimeCharPHC(evttimeminphc,evttimemaxphc,evttimelnphc);
       EarlyActivity(i,earlywph,earlywphpw,earlywphsphere,
                     ph1000,ph500,ph100,lateph1000,lateph500,
                     lateph100,lateph50);
       evtph = ntpEvent->ph.sigcor;

       DupRecoStrips(i,nduprs,dupphrs);

       //figure out batch number
       //note this will only match up with the spillmon info for
       //sure when there are 6 batches
       //if there are only 5 batches, we don't know if there
       //were only 5 batches in the MI, or if the first of those
       //batches went elsewhere.
       batchnumber=FindBatchNumber(evttimemin);


       tree->Fill();
     } // loop over events
   } // loop over entries

   // delete nuparent; // should I do this?
   save.cd();
   pottree->Fill();
   pottree->Write();
   tree->Write();
   save.Write();
   save.Close();

   }

   bool MadTVAnalysis::InFidVol(Int_t evidx){
   if(!LoadEvent(evidx)) return false; //no event found
   return InFidVol();
 }

 bool MadTVAnalysis::InFidVol()
 {
   return InFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,ntpEvent->vtx.z);
 }

 bool MadTVAnalysis::InFidVol(float vx, float vy, float vz){
   // assumes event has been loaded
   // check that vertex is in fiducial volume
   bool is_fid=false;
   // FD: from doc-1351-v2 : there is no coil cut!
   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
         //fiducial criteria on vtx for far detector
     is_fid = true;
     if((vz<0.5 || vz>28.0) ||   //ends
        (vz<16.2 && vz>14.3) ||  //between SMs
        ((vx*vx)+(vy*vy))>14.0) is_fid = false;
   }
   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
     //fiducial criteria on vtx for near detector

     //    float beamzerox = 1.4885;
     //    float beamzeroy = 0.1397-reco_vtxz*TMath::Tan(.058);

     //    float r=sqrt(pow((reco_vtxx-beamzerox),2)+pow((reco_vtxy-beamzeroy),2));


     is_fid = false;
     // MAK: Sept 19, 2005: Modified to account for beam slope
     /*
     double radius = pow(vx-kXcenter,2) 
       + pow(vy-(kYcenter -vz*TMath::Tan(0.058)),2);
     radius=sqrt(radius);
     */

     // MAK: removed beam slope. Dave/Niki don't do it.
     //double radius = sqrt(vx*vx + vy*vy);
     double radius = pow(vx-kXcenter,2) 
       + pow(vy-kYcenter,2);
     radius=sqrt(radius);

     // restrict vertices to target region and 1m radius
     if(vz>=1.0 && vz<=5.0 && radius<=1.0) is_fid = true;    
   }
   return is_fid;
 }
 /*
 Float_t MadTVAnalysis::RecoMKMuEnergy(Int_t& opt,Int_t itrk, bool isdata, const Detector::Detector_t det){
   const float mumass=0.10566; // GeV
   if(LoadTrack(itrk)){
      float mr=ntpTrack->momentum.range;
      mr=EnergyCorrections::CorrectMomentumFromRange(mr,isdata,det);
      float mc=0.0; // signed momentum from curvature
      if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
      mc=EnergyCorrections::CorrectSignedMomentumFromCurvature(mc,isdata,det);

     if(opt==0){  
       //return the most appropriate measure of momentum
       // assign opt based on our choice
       if(ntpTrack->fidall.dr>0.5&&ntpTrack->fidall.dz>0.5) {
         opt=2;
         return sqrt(mr*mr+ mumass*mumass);
       }
       else {
         opt=1;
         // in R1.9 the tracker will apparently return (q/p)=0.0
         // maybe it's when a track looks perfectly rigid?
         // if so, we have to do something
         // I don't want to use the range, that could be very wrong
         // but wrong in a more subtle way ...
         // so, we'll return an obviously ridiculous value of 10 TeV
         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
         else return sqrt(mc*mc+mumass*mumass);
       }
     }
     else if(opt==1) { //return curvature measurement
         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
         else return sqrt(mc*mc + mumass*mumass);
     }
     else if(opt==2) //return range measurement
       return sqrt(mr*mr + mumass*mumass);
     else return 0;
   }
   return 0.;
 }

 Float_t MadTVAnalysis::RecoShwEnergy(int ishw, const Detector::Detector_t& det, int mode, bool isdata){

   Float_t result = 0;
   // Return Jim's calculation
   Bool_t ok = LoadShower(ishw);
   if(!ok) return 0;

   // test for shwph.linCCGeV<=0.0
   // if so, assume that this is an R1.16 object
   // and use ph.gev
   result = ntpShower->shwph.linCCgev;
   if(result<=0.0) result=(ntpShower->ph.gev/1.23);
   else{
      result = EnergyCorrections::CorrectShowerEnergy(result,det,CandShowerHandle::kCC,mode,isdata);
   }
   return result;

 }

Float_t MadTVAnalysis::RecoAnalysisEnergy(Int_t event){
  //  cout<<"A ntpShower: "<<ntpShower<<endl;

  if(!LoadEvent(event)) return 0;
  int largestTrack = -1;
  //  cout<<"B ntpShower: "<<ntpShower<<endl;
  LoadLargestTrackFromEvent(event,largestTrack);
  // cout<<"C ntpShower: "<<ntpShower<<endl;
  int opt = 0;
  float emu = RecoMKMuEnergy(opt,largestTrack);
  // cout<<"D ntpShower: "<<ntpShower<<endl;
  int largestShower = -1;
  LoadLargestShowerFromEvent(event,largestShower);
  float ehad = RecoShwEnergy(largestShower);
  // cout<<"E ntpShower: "<<ntpShower<<endl;
  float ereco = emu + ehad;
  return ereco;

}
*/
Float_t MadTVAnalysis::RecoMuDCosNeuFD(Int_t itr, Float_t* /* vtx */){
  if(!LoadTrack(itr)) return 0.;

  Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
  Float_t bl_y = sqrt(1. - bl_z*bl_z);
  Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
  
  return  costhbl;
}

Float_t MadTVAnalysis::RecoMuDCosNeuND(Int_t itr, Float_t* /* vtx */){
  if(!LoadTrack(itr)) return 0.;
  /*
   // simple correction based on the vertical position
  float vtxy=0;
  if(vtx) vtxy=vtx[1]; // in meters
  // cosine of the typical neutrino angle in the yz plane
  // calculated as py / sqrt( py^2 + pz^2)
  float nu_cos = -5.799E-2;
  if(vtxy>-2.0 && vtxy<2.0){ //prevents further nuttyness if vtxy is silly
    nu_cos += vtxy*1.23304E-3 
      + vtxy*vtxy*1.08212E-5 
      + vtxy*vtxy*vtxy*(-4.634E-5);
  }
  */
  float nu_cos = -5.799E-2;
  float nu_sin = sqrt(1 -nu_cos*nu_cos);
  float cosz = ntpTrack->vtx.dcosz*nu_sin + ntpTrack->vtx.dcosy*nu_cos;

  return cosz;

}

/*
Bool_t MadTVAnalysis::PassBasicCuts(Int_t event){
  return kTRUE;

}
*/

Float_t MadTVAnalysis::ClosestPointToLine(const Float_t& x1, 
                                          const Float_t& y1,
                                          const Float_t& x2,
                                          const Float_t& y2,
                                          const Float_t& xpt,
                                          const Float_t& ypt,
                                          Float_t& x, Float_t& y){
  // given a line with end points (x1,y1) (x2,y2) 
  // and a spatial point (xpt,ypt)
  // report the closest point on the line (x,y) to the spatial point
  // and compute the distance from (xpt,ypt) to (x,y)
  Float_t u= (xpt-x1)*(x2-x1)+(ypt-y1)*(y2-y1);
  u/=sqrt( pow(x2-x1,2) + pow(y2-y1,2));
  x=x1+u*(x2-x1);
  y=y1+u*(y2-y1);
  Float_t dist =sqrt(pow(xpt-x, 2)+pow(ypt-y,2));
  return dist;
}
/*
void MadTVAnalysys::MakeNDPlaneOutline(const PlaneCoverage::PlaneCoverage_t& pc
                                       const PlaneView::PlaneView_t& pv,
                                       vector<Float_t>& v){
  v.clear();
  using PlaneCoverage;
  // I only work on the near detector
  if(!( (pc==kNearPartial) || (pc==kNearFull) || (pc==kTotal) ) ) return;
  using PlaneView;
  // I only work on U and V planes
  if(!( (pv==kU)||(pv==kV)) ) return;
  
  

}
*/

Int_t MadTVAnalysis::NDRadialFidVol(Float_t x,Float_t y, Float_t z){

  Float_t rings[5]={0.1,0.25,0.5,1.0,1.5}; // in meters
  Int_t result=0;
  for (unsigned int i=0; i<5; i++){
    Float_t radius = pow(x-kXcenter,2) 
      + pow(y-(kYcenter - z*TMath::Tan(0.058)),2);
    radius=sqrt(radius);
    // restrict vertices to target region and 1m radius
    if(z>=kVetoEndZ && 
       z<=kTargEndZ &&
       radius<=rings[i]) result|=(1<<i);
  }
  return result;
}

Bool_t MadTVAnalysis::PittInFidVol(Float_t x, Float_t y, Float_t z, Float_t u, Float_t v){

  
  // Is the event vertex, (x,y,u,v,z) inside the fiducial volume
  // used by the pittsburgh group?

  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;

}

Int_t MadTVAnalysis::PittTrkContained(Float_t x, Float_t y, Float_t z, 
                                      Float_t u, Float_t v){
  // what is the containment status of this track 
  // as defined by the pitt group
  //
 // takes the track end point (x,y,u,v,z)
  // returns:
  // 1 = track is fully contained in the upstream region
  // 2 = track is partially contained in the upstream region
  // 3 = track is fully contained in the downstream region
  // 4 = track is partially contained in the downstream region
  // 
  // the rationale is that these categories have different momentum resolution
  //

  const Float_t radsq = x*x + y*y;
  int result=0;
  if(z<7.0) {
    // track in upstream region
    static const Float_t coil_cut=0.8*0.8;
    // was asking for u>3 && u<1.8 here before -- a bug
    if( (u>0.3 && u<1.8) && (v>-1.8 && v<-0.3) && (x<2.4) && (radsq>coil_cut) )
      result=1;
    else result=2;
  }
  else{
    // downstream region
    // uses an ellipse in x,y
    // major axis (x) = a = 1.7 m
    // minor axis (y) = b = 1.4 m
    // centered on x,y = (0.8,0)
    // (x-x0)^2/a^2 + (y-y0)^2/b^2 = 1
    // should change this to 0.8
    //    static const Float_t coil_cut=0.5*0.5; // note differs from above

    //changed from 0.5 to 0.8 by TV on 7-29-05
    static const Float_t coil_cut=0.8*0.8;
    static const Float_t x0=0.8;
    static const Float_t y0=0.0;
    static const Float_t a=1.7;
    static const Float_t b=1.4;
    const Float_t xsc = (x-x0)/a; // rescale ellipse to unit circle
    const Float_t ysc = (y-y0)/b;
    //MAK: cut on z was 16.0 till Aug 2, 2005
    if( (sqrt(xsc*xsc + ysc*ysc)<1.0) && (radsq>coil_cut) && (z<15.6)) 
      result = 3;
    else result = 4;
  }
  return result;
  
}

Int_t MadTVAnalysis::InPartialRegion(UShort_t plane, UShort_t strip){
  // figure out if this (plane,strip) corresponds to something in
  // the partially instrumented region
  //
  // this region is defined as:
  // v planes: (strip<=4 || strip>=67)
  // partial u: (strip==0 || strip=63)
  // full u: (strip<=26 || strip>=88) 
  //
  // if so, return 1
  // if not, return -1
  // if error, return 0


  // make a lookup ptype to hold the type of each plane
  // 1 = v partial   2 = u partial
  // 3 = v full   4 = u full
  // 0 = uninstrumented
  static bool first=true;
  static UShort_t ptype[282];
  if(first){
    ptype[0]=0;
    for(int i=1; i<=281; i++){
      if(i%2==0) ptype[i]=1; // a v plane
      else ptype[i]=2; // a u plane
      if((i-1)%5 == 0) ptype[i]+=2; // fully instrumented
      else if(i>120) ptype[i]=0; // not instrumented
    }
    first=false;
  }
  if(plane>281){
    //    std::cerr<<"InPartialRegion passed plane = "<<plane<<std::endl;
    return 0;
  }
  UShort_t pt = ptype[plane];
  
  Int_t result;
  switch(pt){
  case 1:
  case 3:
    if(strip<=4 || strip>=67) result=1;
    else result = -1;
    break;
  case 2:
    if(strip==0 || strip == 63) result=1;
    else result = -1;
    break;
  case 4:
    if(strip<=26 || strip>=88) result=1;
    else result = -1;
    break;
  case 0:
  default:
    result=0;
    break;
  }
  return result;

}

void MadTVAnalysis::SigInOut(Float_t& sigfull, Float_t& sigpart, 
                             Float_t& trkfull, Float_t& trkpart){
  // compute the amount of signal in the partially instrumented region
  // and the amount in the fully instrumented region
  //
  // this region is defined as:
  // v planes: (strip<=4 || strip>=67)
  // partial u: (strip==0 || strip=63)
  // full u: (strip<=26 || strip>=88) 

  sigfull=sigpart=0;
 
  // loop over all strips in the event
  // and sum the signals in the partial and full regions
  for(int i=0; i<ntpEvent->nstrip; i++){
    if(!LoadStrip(ntpEvent->stp[i])) continue;
    Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip);
    if(pr==1){
      sigpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
    }
    else if(pr==-1){
      sigfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
    }
  }
  
  // loop over those strips in the primary track

  trkfull=trkpart=0;

  if(ntpTrack){
    for(int i=0; i<ntpTrack->nstrip; i++){
      if(!LoadStrip(ntpTrack->stp[i])) continue;
      Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip);
      if(pr==1){
        trkpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
      }
      else if(pr==-1){
        trkfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
      }
    }
  }
  
}

void MadTVAnalysis::SigInOut(Int_t trkidx, Float_t& sigfull, Float_t& sigpart, 
                             Float_t& trkfull, Float_t& trkpart){
  // will try to load the track pointed to by trkidx
  if(trkidx==-1) ntpTrack=0;
  else if (!(LoadTrack(trkidx))) ntpTrack=0;
  SigInOut(sigfull,sigpart,trkfull,trkpart);
}

void MadTVAnalysis::SetBECFile(const char* f){
  if(f){
    fBECConfig.Set("beam:flux_file",f);
  }
}

Int_t MadTVAnalysis::NPlanesInObj(TObject *rec, TObject *obj, float phcut, float& sumph)
{
  sumph=0.0;

  std::set<UShort_t> planes;

  NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec);
  NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec);

  Int_t *stripindex=0;
  Int_t nstrip=0;
  NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj);
  NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj);
  NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj);
  NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj);

  if(slc!=0){
    stripindex=slc->stp;
    nstrip=slc->nstrip;
    //    cout<<"Counting planes in slice"<<endl;
  }
  else if(trk!=0){
    stripindex=trk->stp;
    nstrip=trk->nstrip;
    //    cout<<"Counting planes in trk"<<endl;
  }
  else if(evt!=0){
    stripindex=evt->stp;
    nstrip=evt->nstrip;
    //    cout<<"Counting planes in event"<<endl;
  }
  else if(shw!=0){
    stripindex=shw->stp;
    nstrip=shw->nstrip;
    //    cout<<"Counting planes in event"<<endl;
  }
  
  TClonesArray *striplist=0; 
  int TOTSTRIPS=0;
  if(srrec!=0){
    striplist = srrec->stp;
    TOTSTRIPS=srrec->evthdr.nstrip;
    //    cout<<"Found a sr"<<endl;
  }
  else{
    striplist = strec->stp;
    TOTSTRIPS=strec->evthdr.nstrip;
    //    cout<<"found a st"<<endl;
  }
  
  if(striplist!=0){
    //    cout<<"Looping over strip list, nstrip="<<nstrip<<endl;
    for(int i=0;i<nstrip;i++){
      int sindex = stripindex[i];
      if(sindex<0) continue;
      if(sindex<TOTSTRIPS){
        NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>
          ((*striplist)[sindex]);
        //      cout<<"got strip "<<sindex<<" plane "<<strip->plane
        //          <<" strip "<<strip->strip
        //          <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl;
        if(strip->ph0.pe+strip->ph1.pe>phcut){
          planes.insert(strip->plane);
          sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
        }     
      }
    }
  }
  else{
    //    cout<<"Getting number of planes in snarl"<<endl;
    for(int i=0;i<TOTSTRIPS;i++){
      NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]);
      //      cout<<"got strip "<<i<<" plane "<<strip->plane
      //          <<" strip "<<strip->strip
      //          <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl;
      if(strip->ph0.pe+strip->ph1.pe>phcut){
        planes.insert(strip->plane);
        sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
      }     
    }
  }
    
  //  cout<<"I count "<<planes.size()<<" planes > "<<phcut<<" pe"<<endl;
  return planes.size();
}


Int_t MadTVAnalysis::NStripsInObj(TObject *rec, TObject *obj, float phcut, float& sumph)
{
  sumph=0.0;
  std::vector<UShort_t> strips;

  NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec);
  NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec);

  Int_t *stripindex=0;
  Int_t nstrip=0;
  NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj);
  NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj);
  NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj);
  NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj);

  if(slc!=0){
    stripindex=slc->stp;
    nstrip=slc->nstrip;
  }
  else if(trk!=0){
    stripindex=trk->stp;
    nstrip=trk->nstrip;
  }
  else if(evt!=0){
    stripindex=evt->stp;
    nstrip=evt->nstrip;
  }
  else if(shw!=0){
    stripindex=shw->stp;
    nstrip=shw->nstrip;
  }
  
  TClonesArray *striplist=0; 
  int TOTSTRIPS;
  if(srrec!=0){
    striplist = srrec->stp;
    TOTSTRIPS=srrec->evthdr.nstrip;
  }
  else{
    striplist = strec->stp;
    TOTSTRIPS=strec->evthdr.nstrip;
  }
  
  //  cout<<"CHECKING nstrips="<<nstrip<<endl;
  if(striplist!=0){
    for(int i=0;i<nstrip;i++){
      int sindex=stripindex[i];
      if(sindex<0) continue;
      NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[sindex]);
      if(strip->ph0.pe+strip->ph1.pe>phcut){
        strips.push_back(strip->plane);
        sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
      }     
    }
  }
  else{
    if(srrec!=0){
      nstrip=srrec->evthdr.nstrip;
    }
    else{
      nstrip=strec->evthdr.nstrip;
    }
    for(int i=0;i<nstrip;i++){
      NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]);
      //      cout<<"Checking "<<i<<" stripph: "<<strip->ph0.pe<<" "<<strip->ph1.pe<<endl;
      if(strip->ph0.pe+strip->ph1.pe>phcut){
        strips.push_back(strip->plane);
        sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
      }     
    }
  }
    
  //  cout<<"returning "<<strips.size()<<endl;
  return strips.size();
}

void MadTVAnalysis::GetEvtTimeChar(double &evttimemin, double &evttimemax, double &evttimeln)
{
  //cut and pasted, with minor changes from MadDpAnalysis.cxx
  //  double trgtime=eventSummary->trigtime;
  double trgtime=ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1.e9;
  double timemax=0.;
  double timemin=1.e10;

  //  Find max min time of events by loading strips for ND 
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime=ntpStrip->time1-trgtime;
      if(striptime<=timemin) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
  }      
  
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      //      Double_t striptime1=ntpStrip->time1;
      //      Double_t striptime0=ntpStrip->time0;
      Double_t striptime=(striptime0+striptime1)/2.;
      //      if(striptime1>0 && striptime0<0)  striptime=striptime1;
      //      if(striptime0>0 && striptime1<0)  striptime=striptime0;
      //      if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
      if(striptime1>-1 && striptime0<-1)  striptime=striptime1;
      if(striptime0>-1 && striptime1<-1)  striptime=striptime0;
      if(striptime0>-1 && striptime1>-1)  striptime=(striptime0+striptime1)/2.;

      if(striptime<=timemin) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
  }

    evttimemax=timemax;
    evttimemin=timemin;
    evttimeln=timemax-timemin;

}

void MadTVAnalysis::GetEvtTimeCharPHC(double &evttimemin, double &evttimemax, double &evttimeln)
{
  //same as GetEvtTimeChar but only looking at hits with 
  //pulse height > 200 sig cor (~2pe) in near det
  //pulse height > 160 sig cor in far det
  double trgtime=eventSummary->trigtime;
  double timemax=0.;
  double timemin=1.e10;

  //  Find max min time of events by loading strips for ND 
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime=ntpStrip->time1-trgtime;
      if(ntpStrip->ph1.sigcor>200){
        if(striptime<=timemin) timemin=striptime;
        if(striptime>=timemax) timemax=striptime;
      }
    }
  }      
  
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      Double_t striptime=(striptime0+striptime1)/2.;
      if(striptime1>0 && striptime0<0)  striptime=striptime1;
      if(striptime0>0 && striptime1<0)  striptime=striptime0;
      if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
      if(ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor>160){
        if(striptime<=timemin) timemin=striptime;
        if(striptime>=timemax) timemax=striptime;
      }
    }
  }
  
  evttimemax=timemax;
  evttimemin=timemin;
  evttimeln=timemax-timemin;

}



float MadTVAnalysis::MuonShowerEnergy(const NtpSREvent* evt,const NtpSRTrack* trk){  
  // Purpose: try to sum up brems along the track

  static TH2* hr = new TH2F("h_mushw_rph",
                            "; track - shw r dist (m); shower energy (GeV)", 
                            10,0,0.2,40,0,2);

  static TH2* hz = new TH2F("h_mushw_zph",
                            "; track - shw z dist (m); shower energy (GeV)", 
                            20,-0.2,0.2,40,0,2);

  static TH2* ht = new TH2F("h_mushw_tph",
                            "; track - shw t dist (m); shower energy (GeV)", 
                            600,-1e-6,10e-6,40,0,2);
  

  if(evt==0 || trk==0) return 0;
  float result=0;
  NtpSRShower* save_shwr = ntpShower;  
  // loop through all showers in event
  const float max_vtx_dist=14*0.06; // 2 interaction lengths
  for(int i=0; i<evt->nshower; i++){
    LoadShower(evt->shw[i]);
    if(!ntpShower) continue;
    // calculate distance to vertex
    // don't want to use showers near the vertex
    float dist_vtx = pow(ntpShower->vtx.z - ntpEvent->vtx.z,2)
      + pow(ntpShower->vtx.u - ntpEvent->vtx.u,2)
      + pow(ntpShower->vtx.v - ntpEvent->vtx.v,2);
    dist_vtx=sqrt(dist_vtx);
    if(dist_vtx<max_vtx_dist) continue;
    // loop through track's strip array
    // calculate distance from strip to shower in time and space
    // save spatial and time distance of closest approach to track
    // if shower is close enough add the shower energy to running sum.
    float min_dist=1e6; float min_r=1e6; float min_z=1e6; float min_t=1e6;
    float min_ph=0; float add_ph=0;
    for(int j=0; j<trk->nstrip; j++){
      float dist_trk= pow(trk->stpz[j]-ntpShower->vtx.z,2)
        + pow(trk->stpu[j]-ntpShower->vtx.u,2)
        + pow(trk->stpv[j]-ntpShower->vtx.v,2);
      dist_trk=sqrt(dist_trk);
      float tdist = (trk->vtx.t + (trk->stpz[j]-trk->vtx.z)/3e8) 
                         - ntpShower->vtx.t;
      
        // = trk->stpt0[j]*trk->stpph0[j]+trk->stp1[j]*trk->stpph1[j];
        //      if( (trk->stpph0[j]+trk->stpph1[j]) > 0){
        //      tdist/=(trk->stpph0[j]+trk->stpph1[j]);
        //      }
      
      if(fabs(dist_trk)<fabs(min_dist)) {
        min_ph=ntpShower->shwph.EMgev;
        min_z = fabs(ntpShower->vtx.z-trk->stpz[j]);
        min_r = sqrt( pow(ntpShower->vtx.u-trk->stpu[j],2)
                      + pow(ntpShower->vtx.v-trk->stpv[j],2));
        min_t = tdist;
      }

      const float max_trk_dist=1.76*3.0; // 3 radiation lengths
      const float max_trk_tdist=40e-9;
      if(fabs(dist_trk)<max_trk_dist && fabs(tdist)<max_trk_tdist){
        // shower correllated with track
        add_ph=min_ph;
      }
    }
    result+=min_ph; // 
    
    // fill histograms
    hr->Fill(min_r,min_ph);
    hz->Fill(min_z,min_ph);
    ht->Fill(min_t,min_ph);

  }// end loop over showers

  // reset shower pointer
  ntpShower=save_shwr;
  // all done
  return result;
}

void MadTVAnalysis::EarlyActivity(int evidx, float &earlywph, 
                                  float &earlywphpw, 
                                  float &earlywphsphere, 
                                  float &ph1000, 
                                  float &ph500, 
                                  float &ph100,
                                  float &lateph1000, 
                                  float &lateph500, 
                                  float &lateph100,
                                  float &lateph50)
{
  //zero the variables
  earlywph=0.;
  earlywphpw=0.; 
  earlywphsphere=0.;
  ph1000=0.;
  ph500=0.;
  ph100=0.;
  lateph1000=0.;
  lateph500=0.;
  lateph100=0.;
  lateph50=0.;

  if(!LoadEvent(evidx)) return; //no event found
  
  //find the time of the earliest strip in the event
  double earliestEventTime = 1.e20;
  double latestEventTime = 0.;
  map<int,int> eventPlanes;
  double trigtime = eventSummary->trigtime;

  //loop over strips in the event
  //find the earliest strip time 
  //make a map of the pmt's hit in this event
  for(int s=0;s<ntpEvent->nstrip;s++){
    int stripindex = ntpEvent->stp[s];
    if(!LoadStrip(stripindex)) continue;
    //calculate ph weighted strip time (over the two ends,
    //should be equivalent to ph1 in ND
    float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+
                   ntpStrip->ph0.sigcor*ntpStrip->time0)/
                  (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor));
      
    double time = 1.e9*(phwt-trigtime);
    if(time<earliestEventTime){
      earliestEventTime = time;
    }
    if(time>latestEventTime){
      latestEventTime=time;
    }
    if(ntpStrip->ph1.sigcor>0){
      if(eventPlanes.find(ntpStrip->pmtindex1)==eventPlanes.end()){
        eventPlanes[ntpStrip->pmtindex1] = 1;
      }
    }
    if(ntpStrip->ph0.sigcor>0){
      if(eventPlanes.find(ntpStrip->pmtindex0)==eventPlanes.end()){
        eventPlanes[ntpStrip->pmtindex1] = 1;
      }
    }
  }

  //now find the weigthed sum of ADC for the early strips - make sure the early
  //strips come from the same pmts as the event strips.
  for(unsigned int s=0;s<eventSummary->nstrip;s++){
    if(!LoadStrip(s)) continue;
    //work in ns
    float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+
                   ntpStrip->ph0.sigcor*ntpStrip->time0)/
                  (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor));
    double newtime = 1.e9*(phwt-trigtime);
    double dt=(earliestEventTime-newtime);
    double dtlate=newtime-latestEventTime;
    float d=1000.;
    if((int)ntpStrip->planeview==PlaneView::kU){
      d=sqrt((pow(ntpEvent->vtx.u-ntpStrip->tpos,2)+
              pow(ntpEvent->vtx.z-ntpStrip->z,2)));
    }
    else{
      d=sqrt((pow(ntpEvent->vtx.v-ntpStrip->tpos,2)+
              pow(ntpEvent->vtx.z-ntpStrip->z,2)));
    }
    if(d<fEarlyPlaneSphere&&dtlate>0.){
      if(dtlate<1000.){
        lateph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
      if(dtlate<500.){
        lateph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
      if(dtlate<100.){
        lateph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
      if(dtlate<50.){
        lateph50+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
    }   
    if(d<fEarlyPlaneSphere&&dt>0.){
      if(dt<1000.){
        ph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
      if(dt<500.){
        ph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
      if(dt<100.){
        ph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
      }
    }
    if(dt>0.&&dt<1000.*fEarlyActivityWindowSize){
      double eadc = ((ntpStrip->ph1.sigcor+ntpStrip->ph1.sigcor)*
                     TMath::Exp(-1.*dt/fPMTTimeConstant));
      //if this strip has activity in the same pmt as our event,
      //sum up the early pulse height
      if(eventPlanes.find(ntpStrip->pmtindex1)!=eventPlanes.end()){
        earlywph += eadc;
      }      
      //if this strip is within some radius of event vertex, 
      //add up the early pulse height
      if(d<fEarlyPlaneSphere){
        earlywphsphere+=eadc;
      }
      //if this strip is within the plane window of event vertex,
      //add up the early pulseheight
      if(ntpStrip->plane>=ntpEvent->vtx.plane&&
         ntpStrip->plane<=ntpEvent->vtx.plane+fNPlaneEarlyAct){
        earlywphpw+=eadc;
      }
    }
  }
}


void MadTVAnalysis::FillMCHists(TFile* fout){

  static bool first=true;

  static TH1* h_numu_pf_tot = new TH1F("h_numu_pf_tot","CC; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_pf_qe = new TH1F("h_numu_pf_qe","QE; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_pf_res = new TH1F("h_numu_pf_res","RES; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_pf_dis = new TH1F("h_numu_pf_dis","DIS; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  static TH1* h_numu_pf_nc = new TH1F("h_numu_pf_nc","NC; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  if(first){
    h_numu_pf_tot->SetDirectory(fout);
    h_numu_pf_qe->SetDirectory(fout);
    h_numu_pf_res->SetDirectory(fout);
    h_numu_pf_dis->SetDirectory(fout);
    h_numu_pf_nc->SetDirectory(fout);
  }

  static TH1* h_numubar_pf_tot = new TH1F("h_numubar_pf_tot","CC; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_pf_qe = new TH1F("h_numubar_pf_qe","QE; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_pf_res = new TH1F("h_numubar_pf_res","RES; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_pf_dis = new TH1F("h_numubar_pf_dis","DIS; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  static TH1* h_numubar_pf_nc = new TH1F("h_numubar_pf_nc","NC; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  if(first){
    h_numubar_pf_tot->SetDirectory(fout);
    h_numubar_pf_qe->SetDirectory(fout);
    h_numubar_pf_res->SetDirectory(fout);
    h_numubar_pf_dis->SetDirectory(fout);
    h_numubar_pf_nc->SetDirectory(fout);
  }



  static TH1* h_numu_1m_tot = new TH1F("h_numu_1m_tot","CC; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_1m_qe = new TH1F("h_numu_1m_qe","QE; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_1m_res = new TH1F("h_numu_1m_res","RES; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numu_1m_dis = new TH1F("h_numu_1m_dis","DIS; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  static TH1* h_numu_1m_nc = new TH1F("h_numu_1m_nc","NC; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  if(first){
    h_numu_1m_tot->SetDirectory(fout);
    h_numu_1m_qe->SetDirectory(fout);
    h_numu_1m_res->SetDirectory(fout);
    h_numu_1m_dis->SetDirectory(fout);
    h_numu_1m_nc->SetDirectory(fout);
  }

  static TH1* h_numubar_1m_tot = new TH1F("h_numubar_1m_tot","CC; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_1m_qe = new TH1F("h_numubar_1m_qe","QE; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_1m_res = new TH1F("h_numubar_1m_res","RES; true neutrino energy (GeV)",
                                 240,0.0,120.0);

  static TH1* h_numubar_1m_dis = new TH1F("h_numubar_1m_dis","DIS; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  static TH1* h_numubar_1m_nc = new TH1F("h_numubar_1m_nc","NC; true neutrino energy (GeV)",
                                 240,0.0,120.0);
  if(first){
    h_numubar_1m_tot->SetDirectory(fout);
    h_numubar_1m_qe->SetDirectory(fout);
    h_numubar_1m_res->SetDirectory(fout);
    h_numubar_1m_dis->SetDirectory(fout);
    h_numubar_1m_nc->SetDirectory(fout);
  }


  static TH1* h_other_pf = new TH1F("h_other_pf",
                                    "other; true neutrino energy (GeV)",
                                    240,0.0,120.0);

  static TH1* h_other_1m = new TH1F("h_other_1m",
                                    "other; true neutrino energy (GeV)",
                                    240,0.0,120.0);
  if(first){
    h_other_pf->SetDirectory(fout);
    h_other_1m->SetDirectory(fout);
  }
  if(!mcHeader) {
    first=true;
    return;
  }
  for(Int_t i=0; i<mcHeader->nmc; i++){
    if(!LoadTruth(i)) continue; // no ntpmctruth
    const Float_t k = 1.0/sqrt(2.0);
    Float_t x = ntpTruth->vtxx;
    Float_t y = ntpTruth->vtxy;
    Float_t z = ntpTruth->vtxz;
    Float_t u = k*(y+x);
    Float_t v = k*(y-x);
    Bool_t inpf = PittInFidVol(x,y,z,u,v);    
    Int_t ndr = NDRadialFidVol(x,y,z);
    Bool_t inr =kFALSE;
    if((ndr&0x8)>0) inr=kTRUE;
    Int_t ccnc = ntpTruth->iaction;
    Int_t process = ntpTruth->iresonance;
    Int_t inu = ntpTruth->inu;
    Float_t E = ntpTruth->p4neu[3];

    
    if(inu==14){ // neutrinos
      if(ccnc==0){
        if(inpf) h_numu_pf_nc->Fill(E);
        if(inr) h_numu_1m_nc->Fill(E);
      }
      else if(ccnc==1){
        if(inpf) h_numu_pf_tot->Fill(E);
        if(inr) h_numu_1m_tot->Fill(E);
        switch(process){
        case 1001:
          if(inpf) h_numu_pf_qe->Fill(E);
          if(inr) h_numu_1m_qe->Fill(E);
          break;
        case 1002:
          if(inpf) h_numu_pf_res->Fill(E);
          if(inr) h_numu_1m_res->Fill(E);
          break;
        case 1003:
          if(inpf) h_numu_pf_dis->Fill(E);
          if(inr) h_numu_1m_dis->Fill(E);
          break;
        default:
          break;
        }
      }
    }
    else if(inu==-14){
      if(ccnc==0){
        if(inpf) h_numubar_pf_nc->Fill(E);
        if(inr) h_numubar_1m_nc->Fill(E);
      }
      else if(ccnc==1){
        if(inpf) h_numubar_pf_tot->Fill(E);
        if(inr) h_numubar_1m_tot->Fill(E);
        switch(process){
        case 1001:
          if(inpf) h_numubar_pf_qe->Fill(E);
          if(inr) h_numubar_1m_qe->Fill(E);
          break;
        case 1002:
          if(inpf) h_numubar_pf_res->Fill(E);
          if(inr) h_numubar_1m_res->Fill(E);
          break;
        case 1003:
          if(inpf) h_numubar_pf_dis->Fill(E);
          if(inr) h_numubar_1m_dis->Fill(E);
          break;
        default:
          break;
        }
      }
    }
    else{
      if(inpf) h_other_pf->Fill(E);
      if(inr) h_other_1m->Fill(E);
    }
  }
  
  first=false;
}

int MadTVAnalysis::FarTrkContained(Float_t x, Float_t y, Float_t z)
{
  int is_cont=0;
  //arguments should correspond to track end
  //returns:
  //1 if track is partially contained
  //2 if  track is fully contained
  const Float_t r = sqrt(x*x+y*y);

  if(r<=3.5&&r>.5&&((z>=0.5&&z<=14.5)||(z>=16.5&&z<=29.4))){//contained
    is_cont=2;
  }
  else{//not contained
    is_cont=1;
  }
  return is_cont;
}

float MadTVAnalysis::ComputeLowPHRatio()
{
  float lowphsum=0.;
  float totphsum=0.;

  for(unsigned int i=0;i<GetEventSummary()->nstrip;i++){
    if(!LoadStrip(i)) continue;
    totphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe;
    if(ntpStrip->ph1.pe+ntpStrip->ph0.pe<2){
      lowphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe;
    }
  }

  float result=-1.;
  if(totphsum!=0){
    result = lowphsum/totphsum;
  }

  return result;

}

void MadTVAnalysis::DupRecoStrips(int evidx, int &nduprs, float &dupphrs)
{
  //zero the variables
  nduprs=0;
  dupphrs=0.;

  if(!LoadEvent(evidx)) return; //no event found
  
  map<int,float> sindex;
  for(int i=0;i<ntpEvent->nstrip;i++){
    int si = ntpEvent->stp[i];
    if(!LoadStrip(si)) continue;
    float ph = ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor;
    sindex[ntpEvent->stp[i]]=ph;
  }

  for(int i=0;i<eventSummary->nevent;i++){
    if(i==evidx){
      continue;
    }
    LoadEvent(i);  //load a new event
    for(int j=0;j<ntpEvent->nstrip;j++){
      map<int,float>::iterator it(sindex.find(ntpEvent->stp[j]));
      if(it!=sindex.end()){
        nduprs++;
        dupphrs+=it->second;
      }
    }
  }

  //reload the current event
  LoadEvent(evidx);
}
Int_t MadTVAnalysis::FDRCBoundary(){
  Int_t litag=0;
  Int_t numshower=eventSummary->nshower;
  Float_t allph=eventSummary->ph.raw;
  Int_t plbeg=eventSummary->plane.beg;
  Int_t plend=eventSummary->plane.end;

  if (numshower) {
    if (allph>1e6) litag+=10;

    if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
        ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
        ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
        ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
    if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
        ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
        ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
        ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
  }
  return litag;

}
/*
bool MadTVAnalysis::FDRCBoundary(int p1, int p2){

  // try to implement the readout crate boundary condition dave
  // discusses in doc-1351-v2
  // done as in LIPatternFinderSimple()

  const int kNDetectorSegments = 8;
  const int kMinPlane[kNDetectorSegments] = 
    {1,  65, 129, 193, 250, 314, 378, 442 };
  const int kMaxPlane[kNDetectorSegments] = 
    {64, 128, 192, 248, 313, 377, 441, 485 };
  bool result = false;
  // looks for events in which the start and end plane both correspond to
  // the boundaries of a specific readout crate.
  for(int i=0; i<kNDetectorSegments; i++){
    if((p1==kMinPlane[i] && p2==kMaxPlane[i]) 
       ||(p2==kMinPlane[i] && p1==kMaxPlane[i]) ) result=true;
  }
  return result;
}
*/


Float_t MadTVAnalysis::GetNDCoilCurrent(const VldContext& vc){
  Float_t result = 0.0;
  DbiResultPtr<Dcs_Mag_Near> ptr(vc);
  if(ptr.GetNumRows()){
    const Dcs_Mag_Near* row = ptr.GetRow(0);
    result = row->GetCurrent();
  }
  return result;
}

int MadTVAnalysis::FindBatchNumber(double mintime)
{
  //cuts from peter shannahan
    //Batch 0: 1.7e-6 to 3.0e-6
    //Batch 1: 3.3e-6 to 4.5e-6
    //Batch 2: 4.9e-6 to 6.2e-6
    //Batch 3: 6.5e-6 to 7.8e-6
    //Batch 4: 8.1e-6 to 9.4e-6
    //Batch 5: 9.7e-6 to 11.0e-6

  //mintime already has triggertime subtracted
  double tdiff = mintime*1.e6;
  if(tdiff>1.7&&tdiff<3.0) return 0;
  if(tdiff>3.3&&tdiff<4.5) return 1;
  if(tdiff>4.9&&tdiff<6.2) return 2;
  if(tdiff>6.5&&tdiff<7.8) return 3;
  if(tdiff>8.1&&tdiff<9.4) return 4;
  if(tdiff>9.7&&tdiff<11.0) return 5;
  else return -1;
}

---