#include <set>
#include <vector>
#include "Mad/MadMedQEAnalysis.h"
#include "Mad/NearbyEvents.h"
#include "Mad/Moments.h"
#include "DataUtil/EnergyCorrections.h"
// tricia's beam stuff
//#include "Mad/BeamMonMap.h"
#include "MCReweight/GnumiInterface.h"
#include "MCReweight/NuParent.h"
#include "MCReweight/BeamEnergyCalculator.h"
#include "MCReweight/BeamType.h"
#include "Mad/NewMadQEID.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "SpillTiming/SpillTimeFinder.h"


const Float_t MadMedQEAnalysis::kVetoEndZ = 1.2154;
const Float_t MadMedQEAnalysis::kTargEndZ = 3.5914;
const Float_t MadMedQEAnalysis::kCalorEndZ = 7.1554;
const Float_t MadMedQEAnalysis::kHalfCalorEndZ = kTargEndZ+(kCalorEndZ-kTargEndZ)*0.5;
const Float_t MadMedQEAnalysis::kSpecEndZ = 16.656;
const Float_t MadMedQEAnalysis::kXcenter = 1.4885;
const Float_t MadMedQEAnalysis::kYcenter = 0.1397;

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

//********************************************************
MadMedQEAnalysis::MadMedQEAnalysis(TChain *chainSR,TChain *chainMC,
			     TChain *chainTH,TChain *chainEM)
  : MadAnalysis(chainSR, chainMC, chainTH, chainEM)
{

  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",DetectorType::kNear);
  fBECConfig.Set("beam:low_energy_limit",0.5);
  fBECConfig.Set("beam:high_energy_limit",60.0);

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

}

//********************************************************
MadMedQEAnalysis::MadMedQEAnalysis(JobC *j,string path,int entries) : MadAnalysis(j, path, entries)
{
  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",DetectorType::kNear);
  fBECConfig.Set("beam:low_energy_limit",0.5);
  fBECConfig.Set("beam:high_energy_limit",60.0);

  useRePhFracPdf = true;
  useRecoW2Pdf = true;
  useNshowPdf = true;
  useHpeakPdf = true;
  useNvhsHitsPdf = true;
  useNSubS = false;
  useVhsPh = false;
  useQeEnuKinDif = false;
  useUserDefVariable = false;

}

//********************************************************
MadMedQEAnalysis::~MadMedQEAnalysis()
{

}

//********************************************************
void MadMedQEAnalysis::UseRePhFrac(Bool_t tof){
  useRePhFracPdf = tof;
}

//********************************************************
void MadMedQEAnalysis::UseRecoW2(Bool_t tof){
  useRecoW2Pdf = tof;
}

//********************************************************
void MadMedQEAnalysis::UseNshow(Bool_t tof){
  useNshowPdf = tof;
}

//********************************************************
void MadMedQEAnalysis::UseHpeak(Bool_t tof){
  useHpeakPdf = tof;
}

//********************************************************
void MadMedQEAnalysis::UseNvhsHits(Bool_t tof){
  useNvhsHitsPdf = tof;
}

//********************************************************
void MadMedQEAnalysis::UseNSS(Bool_t tof){
  useNSubS = tof;
}

//********************************************************
void MadMedQEAnalysis::UseVHSPh(Bool_t tof){
  useVhsPh = tof;
}

//********************************************************
void MadMedQEAnalysis::UseEnuKinDif(Bool_t tof){
  useQeEnuKinDif = tof;
}

//********************************************************
void MadMedQEAnalysis::UseUserDefVar(Bool_t tof){
  useUserDefVariable = tof;
}

//********************************************************
void MadMedQEAnalysis::CreatePAN(std::string tag, const char* /*bmonpath*/, const char* QEtag)
{
  
  // is this MC??
  this->GetEntry(0);
  const bool file_is_mc
    =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
  // stuff for bmpt reweighting... from tricia
  // only do this if mc

  GnumiInterface* gnumi=0;  
  if(file_is_mc){
    gnumi=new GnumiInterface();
    
    if(!gnumi->Status()) {
      cout << "Warning: MadMedQEAnalysis::CreatePAN() - Flux files not open."
	   << " Will not be filling NuParent info"
	   << endl;
      cout << "Set environmental variable $GNUMIAUX to point to the "
	   << "directory containing the gnumi files to fix this!"
	   << endl;    
    }
    else{
      const char* gnumidir= getenv("GNUMIAUX");
      cout<<"Read gnumi files from "<<gnumidir<<endl;
    }
  }
  //For Beam Reweighting:
  // must have this here, because it's in ntuple
  NuParent *nuparent = new NuParent();

  
  //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 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
  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_cev;           // fully contained. 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 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_eshw;      // reco shower energy  (shw.ph.gev/1.23)
  Float_t reco_mushw;     // showers along muon track  
  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_dircosneu; // reco track vtx dircos w.r.t neutrino
  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 shower vtx
  Float_t reco_shw_vtxy;      // reco y shower vtx
  Float_t reco_shw_vtxz;      // reco z shower 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 end plane

  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;

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

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

  //  std::cout<<"Second Print"<<std::endl;

  //Trees
  // output file
  std::string savename = tag;
  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("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/s",32000);
  tree->Branch("trgtime",&trgtime,"trgtime/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_cev",&is_cev,"is_cev/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("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_t med_qe_pid;
  Int_t med_qe_pass;
  Float_t medPHfrac;
  //  Float_t avShwPhPP;
  Int_t reHits;
  //Int_t numBigH;
  //Float_t rePH;
  Int_t hPeak;
  Float_t protAngU;
  Float_t protAngV;
  //Float_t evtpdaU;
  //Float_t evtpdaV;
  //Float_t phWithProt;
  //Float_t phAwayProt;
  //Float_t phTransProt;
  Float_t prMomU;
  Float_t prMomV;
  Float_t prMomZ;
  Float_t prE1;
  Float_t prE2; 
  Float_t stdHprMomU;
  Float_t stdHprMomV;
  Float_t stdHprMomZ;
  Float_t stdHprE; 
  Float_t muoMomU;
  Float_t muoMomV;
  Float_t muoMomZ;
  Float_t muoE;
  Float_t stdHmuoMomU;
  Float_t stdHmuoMomV;
  Float_t stdHmuoMomZ;
  Float_t stdHmuoE;
  Float_t muDcosZ;
  Float_t muDcosU;
  Float_t muDcosV;
  Float_t muDcosZerr;
  Float_t muDcosUerr;
  Float_t muDcosVerr;
  //Float_t noPhPdaU;
  //Float_t noPhPdaV;
  //Float_t PHdwNhit;
  //Float_t HTangU;
  //Float_t HTangV;
  Float_t HTshwRmsU;
  Float_t HTshwRmsV;
  //Float_t HTvtxDfU;
  //Float_t HTvtxDfV;
  //Int_t medNumShw;
  Float_t stdHprAngU;
  Float_t stdHprAngV;
  Float_t VhsAngU;
  Float_t VhsAngV;
  Float_t allHangU;
  Float_t allHangV;

  Float_t evTotPh;
  Float_t evRePh;
  Float_t phIn0_5;
  Float_t phIn1_0;
  Float_t phIn1_5;
  Float_t phIn2_0;
  Int_t nhIn0_5;
  Int_t nhIn1_0;
  Int_t nhIn1_5;
  Int_t nhIn2_0;

  Float_t teIn0_5;
  Float_t teIn1_0;
  Float_t teIn1_5;
  Float_t teIn2_0;

  Float_t totalTe;

  Float_t userVarble;

  Float_t phInPrange;
  Float_t angSpreadShw;

  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);
  tree->Branch("medPHfrac",&medPHfrac,"medPHfrac/F",32000);
  //tree->Branch("avShwPhPP",&avShwPhPP,"avShwPhPP/F",32000);
  tree->Branch("reHits",&reHits,"reHits/I",32000);
  //tree->Branch("numBigH",&numBigH,"numBigH/I",32000);
  //tree->Branch("rePH",&rePH,"rePH/F",32000);
  tree->Branch("hPeak",&hPeak,"hPeak/I",32000);
  tree->Branch("protAngU",&protAngU,"protAngU/F",32000);
  tree->Branch("protAngV",&protAngV,"protAngV/F",32000);
  //tree->Branch("evtpdaU",&evtpdaU,"evtpdaU/F",32000);
  //tree->Branch("evtpdaV",&evtpdaV,"evtpdaV/F",32000);
  //tree->Branch("phWithProt",&phWithProt,"phWithProt/F",32000);
  //tree->Branch("phAwayProt",&phAwayProt,"phAwayProt/F",32000);
  //tree->Branch("phTransProt",&phTransProt,"phTransProt/F",32000);
  tree->Branch("prMomU",&prMomU,"prMomU/F",32000);
  tree->Branch("prMomV",&prMomV,"prMomV/F",32000);
  tree->Branch("prMomZ",&prMomZ,"prMomZ/F",32000);
  tree->Branch("prE1",&prE1,"prE1/F",32000);
  tree->Branch("prE2",&prE2,"prE2/F",32000);
  tree->Branch("stdHprMomU",&stdHprMomU,"stdHprMomU/F",32000);
  tree->Branch("stdHprMomV",&stdHprMomV,"stdHprMomV/F",32000);
  tree->Branch("stdHprMomZ",&stdHprMomZ,"stdHprMomZ/F",32000);
  tree->Branch("stdHprE",&stdHprE,"stdHprE/F",32000);
  tree->Branch("muoMomU",&muoMomU,"muoMomU/F",32000);
  tree->Branch("muoMomV",&muoMomV,"muoMomV/F",32000);
  tree->Branch("muoMomZ",&muoMomZ,"muoMomZ/F",32000);
  tree->Branch("muoE",&muoE,"muoE/F",32000);
  tree->Branch("stdHmuoMomU",&stdHmuoMomU,"stdHmuoMomU/F",32000);
  tree->Branch("stdHmuoMomV",&stdHmuoMomV,"stdHmuoMomV/F",32000);
  tree->Branch("stdHmuoMomZ",&stdHmuoMomZ,"stdHmuoMomZ/F",32000);
  tree->Branch("stdHmuoE",&stdHmuoE,"stdHmuoE/F",32000);
  tree->Branch("muDcosZ",&muDcosZ,"muDcosZ/F",32000);
  tree->Branch("muDcosU",&muDcosU,"muDcosU/F",32000);
  tree->Branch("muDcosV",&muDcosV,"muDcosV/F",32000);
  tree->Branch("muDcosZerr",&muDcosZerr,"muDcosZerr/F",32000);
  tree->Branch("muDcosUerr",&muDcosUerr,"muDcosUerr/F",32000);
  tree->Branch("muDcosVerr",&muDcosVerr,"muDcosVerr/F",32000);
  //tree->Branch("noPhPdaU",&noPhPdaU,"noPhPdaU/F",32000);
  //tree->Branch("noPhPdaV",&noPhPdaV,"noPhPdaV/F",32000);
  //tree->Branch("PHdwNhit",&PHdwNhit,"PHdwNhit/F",32000);
  //tree->Branch("HTangU",&HTangU,"HTangU/F",32000);
  //tree->Branch("HTangV",&HTangV,"HTangV/F",32000);
  tree->Branch("HTshwRmsU",&HTshwRmsU,"HTshwRmsU/F",32000);
  tree->Branch("HTshwRmsV",&HTshwRmsV,"HTshwRmsV/F",32000);
  //tree->Branch("HTvtxDfU",&HTvtxDfU,"HTvtxDfU/F",32000);
  //tree->Branch("HTvtxDfV",&HTvtxDfV,"HTvtxDfV/F",32000);
  //tree->Branch("medNumShw",&medNumShw,"medNumShw/I",32000);

  tree->Branch("stdHprAngU",&stdHprAngU,"stdHprAngU/F",32000);
  tree->Branch("stdHprAngV",&stdHprAngV,"stdHprAngV/F",32000);
  tree->Branch("VhsAngU",&VhsAngU,"VhsAngU/F",32000);
  tree->Branch("VhsAngV",&VhsAngV,"VhsAngV/F",32000);
  tree->Branch("allHangU",&allHangU,"allHangU/F",32000);
  tree->Branch("allHangV",&allHangV,"allHangV/F",32000);

  tree->Branch("evTotPh",&evTotPh,"evTotPh/F",32000);
  tree->Branch("evRePh",&evRePh,"evRePh/F",32000);
  tree->Branch("phIn0_5",&phIn0_5,"phIn0_5/F",32000);
  tree->Branch("phIn1_0",&phIn1_0,"phIn1_0/F",32000);
  tree->Branch("phIn1_5",&phIn1_5,"phIn1_5/F",32000);
  tree->Branch("phIn2_0",&phIn2_0,"phIn2_0/F",32000);
  tree->Branch("nhIn0_5",&nhIn0_5,"nhIn0_5/I",32000);
  tree->Branch("nhIn1_0",&nhIn1_0,"nhIn1_0/I",32000);
  tree->Branch("nhIn1_5",&nhIn1_5,"nhIn1_5/I",32000);
  tree->Branch("nhIn2_0",&nhIn2_0,"nhIn2_0/I",32000);

  tree->Branch("teIn0_5",&teIn0_5,"teIn0_5/F",32000);
  tree->Branch("teIn1_0",&teIn1_0,"teIn1_0/F",32000);
  tree->Branch("teIn1_5",&teIn1_5,"teIn1_5/F",32000);
  tree->Branch("teIn2_0",&teIn2_0,"teIn2_0/F",32000);

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

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

  tree->Branch("phInPrange",&phInPrange,"phInPrange/F",32000);
  tree->Branch("angSpreadShw",&angSpreadShw,"angSpreadShw/F",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);
 
  Float_t shw_pur, shw_comp;
  tree->Branch("shw_pur",&shw_pur,"shw_pur/F",32000);
  tree->Branch("shw_comp",&shw_comp,"shw_comp/F",32000);

  tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/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("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("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("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");

  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");

  //demux status
  tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b");
  tree->Branch("lowphrat",&lowphrat,"lowphrat/F");
  
  // beam monitoring variables
  Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill;
  //goodbeam==1 if spill meets criteria defined below 
  //(see the line where goodbeam gets set)
  //beamcnf is an int as defined in BeamMonSpill::StatusBits
  UInt_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;

  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("hornI",&hornI,"hornI/D");
  tree->Branch("closestspill",&closestspill,"closestspill/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");

  //coil status
  Short_t coil_stat;
  tree->Branch("coil_stat",&coil_stat,"coil_stat/S");

  //for bmpt reweighting
  tree->Branch("NuParent","NuParent",&nuparent,8000,1);

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

  //  map<VldTimeStamp, BeamMonTV>* beam_mon_map=0;
  //  if(!file_is_mc) beam_mon_map = new BeamMonMap::MakeBeamMonMap(bmonpath);
  //map<VldTimeStamp, BeamMonTV> beam_mon_map
  //=BeamMonMap::MakeBeamMonMap(bmonpath);
  //  std::cout<<"map file: "<<bmonpath<<std::endl;
  //  std::cout<<"map size: "<<beam_mon_map.size()<<std::endl;

  ////////////// beam energy reweighting //////////////////////////////////
  ///// Note: to present a consistent interface the beam energy calculator
  ///// employs a registry.  On an event by event basis, however,
  ///// it directly inputs the configuration paramters to the weighting
  ///// routine
  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
  NewMadQEID qeid; 
  Bool_t fal = false;
  Bool_t tru = true;
  if(useRePhFracPdf) qeid.UseRePHfrac(tru);
  if(!useRePhFracPdf) qeid.UseRePHfrac(fal);
  if(useRecoW2Pdf) qeid.UseRecoWsqd(tru);
  if(!useRecoW2Pdf) qeid.UseRecoWsqd(fal);
  if(useNshowPdf) qeid.UseNshows(tru);
  if(!useNshowPdf) qeid.UseNshows(fal);
  if(useHpeakPdf) qeid.UseHoughPeak(tru);
  if(!useHpeakPdf) qeid.UseHoughPeak(fal);
  if(useNvhsHitsPdf) qeid.UseNVHShits(tru);
  if(!useNvhsHitsPdf) qeid.UseNVHShits(fal);
  if(useNSubS) qeid.UseNsubshows(tru);
  if(!useNSubS) qeid.UseNsubshows(fal);
  if(useVhsPh) qeid.UseVHSPH(tru);
  if(!useVhsPh) qeid.UseVHSPH(fal);
  if(useQeEnuKinDif) qeid.UseQEkinEnuDif(tru);
  if(!useQeEnuKinDif) qeid.UseQEkinEnuDif(fal);
  if(useUserDefVariable) qeid.UseUserVar(tru);
  if(!useUserDefVariable) qeid.UseUserVar(fal);

  qeid.ReadPDFs(QEtag);

  ///////////////////////////// main loop over snarls ////////////////////
    //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);


    //if(!checkeddb){
      //      ntpHeader->GetVldContext().Print();
      //      cout<<"IM HERE!!!!"<<endl;
      //DbiResultPtr<BeamMonSpill> testbs(ntpHeader->GetVldContext(), 
    //Dbi::kDefaultTask, 
    //					Dbi::kDisabled);
      //      cout<<"tried to get the dbiresult ptr"<<endl;
      //checkeddb=true;
      //      cout<<"****************checking for beam monitoring db:"<<endl;
      //if(testbs.GetNumRows()>0){
  //beamdb=true;
  //cout<<"Found the Beam Monitoring Database Tables"<<endl;
  //}
  //  else{	
  //cout<<"Didn't find the Beam Monitoring Database Tables"<<endl;
  //cout<<"Resorting to old beamsummary ntuples"<<endl;
  //  }
  //}

    nevent = eventSummary->nevent;    
    run = ntpHeader->GetRun();
    if(run!=lastrun){
      NRUNS++;
      lastrun = run;
    }
    snarl = ntpHeader->GetSnarl();
    subrun = ntpHeader->GetSubRun();
    trgsrc = ntpHeader->GetTrigSrc();
    spilltype = ntpHeader->GetRemoteSpillType();
    trgtime = eventSummary->trigtime;
    litime = eventSummary->litime;
    coil_stat = detStatus->coilstatus;
    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();

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

    detector = -1; 
    //get detector type:
    const DetectorType::Detector_t det 
      = ntpHeader->GetVldContext().GetDetector();
    if(det==DetectorType::kFar){
      detector = 2;
    }
    else if(det==DetectorType::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
    //////////// BEAM MONITORING ////////////////////////////////
    hbw=vbw=hpos1=vpos1=hpos2=vpos2=hornI=closestspill=0.0;
    beamcnf=0;
    nuTarZ=0.0;
    goodbeam=0;
    horn_on=target_in=n_batches=0;
    tor101=tr101d=tortgt=trtgtd=0.0;
    hadmon=mumon1=mumon2=mumon3=0.0;  

    if(!file_is_mc){
      VldContext vc = ntpHeader->GetVldContext();
      VldTimeStamp vts = vc.GetTimeStamp();
      //      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]);
	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;
      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;
      //      }
      /*
      else{
	// could the following return a reference?
	BeamMonTV btv = BeamMonMap::FindClosestSpill(beam_mon_map,vts);
	hbw=btv.hbw;
	vbw=btv.vbw;
	hpos2=btv.hpos2;
	vpos2=btv.vpos2;
	hornI=btv.hornI;
	closestspill=btv.closestspill;
	
	nuTarZ=btv.nuTarZ;
	//approximately assign beamcnf based on nuTarZ
	if(nuTarZ<20000&&nuTarZ>=0){
	  beamcnf=1;
	  target_in=1;
	}
	else if(nuTarZ>20000&&nuTarZ<60000){
	  beamcnf=4;
	  target_in=1;
	}
	else{
	  beamcnf=5;
	  target_in=1;
	}
	if(fabs(hornI)>10){
	  horn_on=1;
	}
	n_batches=0;
	tor101=0.;
	tr101d=0.;	  
	tortgt=btv.bI;
	trtgtd=0.;
	hadmon=0.;
	mumon1=0.;
	mumon2=0.;
	mumon3=0.;
      }
      */
      //set goodbeam==1 if it passes some basic beam monitoring cuts
      //changing widths to deal with rms instead of fits
      //getting rid of horn_on and target_in since they sometime
      //go haywire
      //	if(horn_on&&target_in&&tortgt>0.1&&hbw<1.5&&vbw<2&&
      //	   hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(closestspill)<2){
      /*
      if(tor101>0.1&&tortgt>0.1&&hbw<2.9&&vbw<2.9&&
	 hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(hornI)>50&&
	 fabs(closestspill)<2){
	goodbeam=1;
      }
      */
      //use Mark D. BMSpillAna to set good beam
      bmana.SetSpill(*bs);
      bmana.SetTimeDiff(closestspill);
      if(bmana.SelectSpill()){
	goodbeam=1;
      }
      else{
	goodbeam=0;
      }

      // attempt to deduce beam type from beam monitoring
      current_beam = BeamType::FromBeamMon(beamcnf);
      
      if(goodbeam==1&&spilltype!=3){//don't count fake spills
	pot+=tortgt;
	pottor101+=tor101;
	pottortgt+=tortgt;
      }
    }

    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

      //cout << "Filling for event: " << i << endl;

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

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

      mcevent = -1;
      is_fid = 0; is_cev = 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_vtx_eshw=0; reco_vtx_wtd_eshw=0;
      reco_qe_enu = 0;   reco_mushw=0;
      reco_qe_q2 = 0; reco_dircosneu = 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;
      // kinematic quantities always positive
      // convention: -1 indicates unfilled variable (NC events)
      reco_y = reco_y = reco_Q2 = reco_W2 = -1;
      
      shw_pur = shw_comp = -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;
      ntotss=0;
      for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;}

      niki_cc_pid=-999;
      dave_cc_pid=-999;

      med_qe_pid=-1001.0;
      med_qe_pass=-1;
      medPHfrac=-2.0;
      //avShwPhPP=-2.0;
      reHits=-2;
      //numBigH=-2;
      //rePH=-2.0;
      hPeak=-2;
      protAngU=-12.0;
      protAngV=-12.0;
      //evtpdaU=-3.0;
      //evtpdaV=-3.0;
      //phWithProt=-2.0;
      //phAwayProt=-2.0;
      //phTransProt=-2.0;
      prMomU=-101.0;
      prMomV=-101.0;
      prMomZ=-101.0;
      prE1=-2.0;
      prE2=-2.0;
      stdHprMomU=-101.0;
      stdHprMomV=-101.0;
      stdHprMomZ=-101.0;
      stdHprE=-2.0;
      muoMomU=-101.0;
      muoMomV=-101.0;
      muoMomZ=-101.0;
      muoE=-2.0;
      stdHmuoMomU=-101.0;
      stdHmuoMomV=-101.0;
      stdHmuoMomZ=-101.0;
      stdHmuoE=-2.0;
      muDcosZ = -100.0;
      muDcosU = -100.0;
      muDcosV = -100.0;
      muDcosZerr = -100.0;
      muDcosUerr = -100.0;
      muDcosVerr = -100.0;
      //noPhPdaU=-3.0;
      //noPhPdaV=-3.0;
      //PHdwNhit=-2.0;
      //HTangU=-12.0;
      //HTangV=-12.0;
      HTshwRmsU=-12.0;
      HTshwRmsV=-12.0;
      //HTvtxDfU=-12.0;
      //HTvtxDfV=-12.0;
      //medNumShw=-2;

      stdHprAngU = -12.0;
      stdHprAngV = -12.0;
      VhsAngU = -12.0;
      VhsAngV = -12.0;
      allHangU = -12.0;
      allHangV = -12.0;


      evTotPh = -2.0;
      evRePh = -2.0;
      phIn0_5 = -2.0;
      phIn1_0 = -2.0;
      phIn1_5 = -2.0;
      phIn2_0 = -2.0;
      nhIn0_5 = -2;
      nhIn1_0 = -2;
      nhIn1_5 = -2;
      nhIn2_0 = -2;

      teIn0_5 = -2.0;
      teIn1_0 = -2.0;
      teIn1_5 = -2.0;
      teIn2_0 = -2.0;

      totalTe = -2.0;

      userVarble = -2.0;

      phInPrange = -2.0;
      angSpreadShw = -12.0;

      //////////////////////////////////////////////////////////////////////


      is_fid = InFidVol();

      //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);
	/*
	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;
      }

      if(PassBasicCuts(event)) {
	//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);

	pass_track = PassTrackCuts(track_index);
	//methods should all be protected against index = -1
	
	reco_emu          = RecoMKMuEnergy(emu_meth,track_index);

	int shower_index  = -1;

	// load the largest shower
	/*
	if(track_index!=-1) {
	  LoadShowerAtTrackVertex(i,track_index,shower_index);
	}
	// shower with the best match to the vertex	
	else LoadLargestShowerFromEvent(i,shower_index);
	*/
	/*
	if(nshower>0){
	  shower_index=ntpEvent->shw[0];
	  LoadShower(shower_index);
	}
	*/

	// use the jim shower!
	LoadShower_Jim(i,shower_index,det);
	reco_eshw         = RecoShwEnergy(shower_index, det,1,!file_is_mc);
	if(shower_index>-1) reco_wtd_eshw = ntpShower->shwph.wtCCgev;
	reco_vtx_eshw = reco_eshw;
	reco_vtx_wtd_eshw = reco_wtd_eshw;

	// figure shower truth business out
	if(is_mc && shower_index>-1){
	  LoadTHShower(shower_index);
	  assert(ntpTHShower);
	  shw_pur=ntpTHShower->purity;
	  shw_comp=ntpTHShower->completeslc;
	}

	// 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 is only useful for FD events
	  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
	  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;
	    trkqp             = ntpTrack->momentum.qp;
	    if(trkqp!=0){
	      trkqp = 1.0/CorrectSignedMomentumFromCurvature(1.0/trkqp,
							     !file_is_mc, det);
	    }

	    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;
	    // 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
	    if(detector==1){//if it's the near detector
	      if(pitt_evt_class>0){ // if it was classified
		int do_meth=0; // emu reco method we should use
		if( (pitt_evt_class == 1) || (pitt_evt_class == 3) ) 
		  do_meth=2; // stoppers --> range
		else if( (pitt_evt_class == 2) || (pitt_evt_class == 4) ) 
		  do_meth=1; // punch-through --> curvature
		if(do_meth==1 || do_meth==2){
		   reco_emu = RecoMKMuEnergy(do_meth,track_index,
					     !file_is_mc,det);
		   reco_enu = reco_emu+reco_eshw;
		}
		emu_meth=do_meth; // was forgetting about this
	      }
	    }
	    else if(detector==2){//if it's the near detector
	      int do_meth=FarTrkContained(ntpTrack->end.x,
					  ntpTrack->end.y,
					  ntpTrack->end.z);
	      if(do_meth==1 || do_meth==2){
		reco_emu = RecoMKMuEnergy(do_meth,track_index);
		reco_enu = reco_emu+reco_eshw;
	      }
	      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
	    if( (ntpTrack->fit.pass>0) 
		&& (ntpTrack->fit.chi2/ntpTrack->fit.ndof <20 )
		&& abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
	      pitt_trk_qual=1;
	    
	    ////////// parton model/DIS kinematics quantities///////// 
	    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);
	    ///////////////// END kinematics /////////////////////////

	    ///////// recalculate QE energy and Q2 variables /////////
	    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);
	    //////////////// END QE energy Q2 ////////////////////////

	    //marks qeid
	    //med_qe_pid=qeid.CalcQEID(this,i,reco_enu,reco_qe_enu,reco_emu,reco_W2);
	    //med_qe_pass=qeid.PassMedQECut(reco_enu,med_qe_pid);

	    if(qeid.CalcQEVars(this,i,reco_enu,reco_emu)){
	    }

	    medPHfrac=qeid.GetRePHfrac();
	    //avShwPhPP=qeid.GetAvShwPHpPlane();
	    //reHits=qeid.GetReNumHits();
	    reHits=qeid.GetNVHShits();
	    //numBigH=qeid.GetNumBigHits();
	    //rePH=qeid.GetPhRem();
	    hPeak=qeid.GetHoughPeak();
	    protAngU=qeid.GetProtonAngU();
	    protAngV=qeid.GetProtonAngV();
	    //evtpdaU=qeid.GetEvtPdaU();
	    //evtpdaV=qeid.GetEvtPdaV();
	    //phWithProt=qeid.GetPhProjWithProt();
	    //phAwayProt=qeid.GetPhProjAwayProt();
	    //phTransProt=qeid.GetPhProjTransProt();
	    /*
	    prMomU=qeid.GetProtnMomU();
	    prMomV=qeid.GetProtnMomV();
	    prMomZ=qeid.GetProtnMomZ();
	    prE=qeid.GetProtnE();
	    stdHprMomU=qeid.GetStdHepPmomU();
	    stdHprMomV=qeid.GetStdHepPmomV();
	    stdHprMomZ=qeid.GetStdHepPmomZ();
	    stdHprE=qeid.GetStdHepPE();
	    muoMomU=qeid.GetMuMomU();
	    muoMomV=qeid.GetMuMomV();
	    muoMomZ=qeid.GetMuMomZ();
	    muoE=qeid.GetMuE();
	    stdHmuoMomU=qeid.GetStdHMuMomU();
	    stdHmuoMomV=qeid.GetStdHMuMomV();
	    stdHmuoMomZ=qeid.GetStdHMuMomZ();
	    stdHmuoE=qeid.GetStdHMuE();
	    */

	    prMomU=qeid.GetProtonUcomp();
	    prMomV=qeid.GetProtonVcomp();
	    prMomZ=qeid.GetProtonZcomp();
	    prE1=qeid.GetProtonEfromEvent();
	    prE2=qeid.GetProtonEfromMassMom();
	    stdHprMomU=qeid.GetStdHepProtUcomp();
	    stdHprMomV=qeid.GetStdHepProtVcomp();
	    stdHprMomZ=qeid.GetStdHepProtZcomp();
	    stdHprE=qeid.GetStdHepProtE();
	    muoMomU=qeid.GetMuonUcomp();
	    muoMomV=qeid.GetMuonVcomp();
	    muoMomZ=qeid.GetMuonZcomp();
	    muoE=qeid.GetMuonE();
	    stdHmuoMomU=qeid.GetStdHepMuUcomp();
	    stdHmuoMomV=qeid.GetStdHepMuVcomp();
	    stdHmuoMomZ=qeid.GetStdHepMuZcomp();
	    stdHmuoE=qeid.GetStdHepMuE();

	    muDcosZ=ntpTrack->vtx.dcosz;
	    muDcosU=ntpTrack->vtx.dcosu;
	    muDcosV=ntpTrack->vtx.dcosv;
	    muDcosZerr=ntpTrack->vtx.edcosz;
	    muDcosUerr=ntpTrack->vtx.edcosu;
	    muDcosVerr=ntpTrack->vtx.edcosv;
	    //noPhPdaU=qeid.GetNoPHpdaU();
	    //noPhPdaV=qeid.GetNoPHpdaV();
	    //PHdwNhit=qeid.GetPHdwNhits();
	    //HTangU=qeid.GetHoughAngU();
	    //HTangV=qeid.GetHoughAngV();
	    //HTshwRmsU=qeid.GetShwRms75U();
	    //HTshwRmsV=qeid.GetShwRms75V();
	    //HTvtxDfU=qeid.GetHTvtxDifU();
	    //HTvtxDfV=qeid.GetHTvtxDifV();
	    //medNumShw=qeid.GetNshows();

	    HTshwRmsU=qeid.GetHoughRmsU();
	    HTshwRmsV=qeid.GetHoughRmsV();

	    stdHprAngU=qeid.GetStdHepProtAngU();
	    stdHprAngV=qeid.GetStdHepProtAngV();
	    VhsAngU=qeid.GetVHSangU();
	    VhsAngV=qeid.GetVHSangV();
	    allHangU=qeid.GetAllHitAngU();
	    allHangV=qeid.GetAllHitAngV();

	    evTotPh=qeid.GetEvTotPH();
	    evRePh=qeid.GetEvRePH();
	    phIn0_5=qeid.GetPhIn0_5sigCone();
	    phIn1_0=qeid.GetPhIn1_0sigCone();
	    phIn1_5=qeid.GetPhIn1_5sigCone();
	    phIn2_0=qeid.GetPhIn2_0sigCone();
	    nhIn0_5=qeid.GetNhIn0_5sigCone();
	    nhIn1_0=qeid.GetNhIn1_0sigCone();
	    nhIn1_5=qeid.GetNhIn1_5sigCone();
	    nhIn2_0=qeid.GetNhIn2_0sigCone();

	    teIn0_5=qeid.GetTeIn0_5sigCone();
	    teIn1_0=qeid.GetTeIn1_0sigCone();
	    teIn1_5=qeid.GetTeIn1_5sigCone();
	    teIn2_0=qeid.GetTeIn2_0sigCone();

	    totalTe=qeid.GetTeTot();

	    userVarble=qeid.GetUserVariable();

            phInPrange=qeid.GetPhInProtRange();
            angSpreadShw=qeid.GetShwAngSpread();

	    

	    ////////////////// niki's cc 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	    
	    ////////////////// END niki's cc pid /////////////////////////

	    ///////////////////// dave's cc pid //////////////////////////
	    if(dpid.ChoosePDFs(det,current_beam)){
	      dave_cc_pid = dpid.CalcPID(this,event,0);
	    }
	    //////////////////// END dave cc pid /////////////////////////

	    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;


	    ////////////////// access subshowers////////////
	    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
	    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();
	    
	    ///////// done with subshowers ////////////////
	  }

	  

	  if(ntrack==1 && reco_enu>0.0 && reco_enu<140
	     && is_pitt_fid && pitt_trk_qual) pass=1;

	  //if(process==1001 && cc_nc==1 && pass==1){
	  //  cout << "QE event,  PHfrac = " << medPHfrac << endl;
	  //  if(medPHfrac==-1){
	  //    cout << "   remaining PH  :" << evRePh << endl;
	  //    cout << "   total event PH:" << evTotPh << endl;
	  //  }
	    //cout << "----------------------" << endl;
	    //cout << "Run:          " << run << endl;
	    //cout << "Snarl:        " << snarl << endl;
	    //cout << "CC or NC:     " << cc_nc << endl;
	    //cout << "Process:      " << process << endl;
	    //cout << "E_neutrino:   " << reco_enu << endl;
	    //cout << "E_muon:       " << reco_emu << endl;
	    //cout << endl;
	    //cout << "p_mom_U:      " << prMomU << endl;
	    //cout << "  truth:      " << stdHprMomU << endl;
	    //cout << "p_mom_V:      " << prMomV << endl;
	    //cout << "  truth:      " << stdHprMomV << endl;
	    //cout << "p_mom_Z:      " << prMomZ << endl;
	    //cout << "  truth:      " << stdHprMomZ << endl;
	    //cout << "p_energy:     " << prE << endl;
	    //cout << "   truth:     " << stdHprE << endl;
	  //}
	  //if(process==1002 && cc_nc==1 && pass==1){
	  //  cout << "RES event, PHfrac = " << medPHfrac << endl;
	  //  if(medPHfrac==-1){
	  //    cout << "   remaining PH  :" << evRePh << endl;
	  //    cout << "   total event PH:" << evTotPh << endl;
	  //  }
	  // }
	  //if(process==1003 && cc_nc==1 && pass==1){
	  //  cout << "DIS event, PHfrac = " << medPHfrac << endl;
	  //  if(medPHfrac==-1){
	  //    cout << "   remaining PH  :" << evRePh << endl;
	  //    cout << "   total event PH:" << evTotPh << endl;
	  //  }
	  // }
	  //if(cc_nc==0 && pass==1){
	  //  cout << "NC event,  PHfrac = " << medPHfrac << endl;
	  //  if(medPHfrac==-1){
	  //    cout << "   remaining PH  :" << evRePh << endl;
	  //    cout << "   total event PH:" << evTotPh << endl;
	  //  }
	  //}
	  
	  

	  //marks qeid
	  med_qe_pid=qeid.CalcQEID(this,i,reco_enu,reco_qe_enu,reco_emu,reco_W2,ntotss);
	  med_qe_pass=qeid.PassMedQECut(reco_enu,med_qe_pid);

	} // 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);

      ///////////////////////// BMPT reweighting /////////////////////////////
      if(file_is_mc){
	if(gnumi->Status()) {
	  if(isST) gnumi->GetParent(strecord,mcevent,*nuparent);
	  else gnumi->GetParent(mcrecord,mcevent,*nuparent);
	}		
      }
      
      //if(pass==1) cout << "Have got variables and pass=1, eg: medPHfrac = " << medPHfrac << endl;

      tree->Fill();

      //if(pass==1) cout << "Filled tree..." << endl;

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

}


//********************************************************
void MadMedQEAnalysis::CreateQEPDFs(const char* QEtag)
{

  // is this MC??
  this->GetEntry(0);
  const bool file_is_mc
    =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
  
  if(!file_is_mc) return;

  //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
  
  //Reco Quantities
  Int_t detector;         // Near = 1, Far = 2;
  Int_t run;              // run number
  Int_t snarl;            // snarl number
 
  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_cev;           // fully contained. 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)

  Int_t pass_track;       // the primary track passes track cuts
  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_eshw;      // reco shower energy  (shw.ph.gev/1.23)
  Float_t reco_mushw;     // showers along muon track  
  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_dircosneu; // reco track vtx dircos w.r.t neutrino
  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 shower vtx
  Float_t reco_shw_vtxy;      // reco y shower vtx
  Float_t reco_shw_vtxz;      // reco z shower 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 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

  // 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 end plane

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

  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;

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

  Float_t etu, etv, elu, elv;
  Float_t swu,swv,wswu,wswv; 

  // mark's QE id
  NewMadQEID qeid; 
  Bool_t fal = false;
  Bool_t tru = true;
  if(useRePhFracPdf) qeid.UseRePHfrac(tru);
  if(!useRePhFracPdf) qeid.UseRePHfrac(fal);
  if(useRecoW2Pdf) qeid.UseRecoWsqd(tru);
  if(!useRecoW2Pdf) qeid.UseRecoWsqd(fal);
  if(useNshowPdf) qeid.UseNshows(tru);
  if(!useNshowPdf) qeid.UseNshows(fal);
  if(useHpeakPdf) qeid.UseHoughPeak(tru);
  if(!useHpeakPdf) qeid.UseHoughPeak(fal);
  if(useNvhsHitsPdf) qeid.UseNVHShits(tru);
  if(!useNvhsHitsPdf) qeid.UseNVHShits(fal);
  if(useNSubS) qeid.UseNsubshows(tru);
  if(!useNSubS) qeid.UseNsubshows(fal);
  if(useVhsPh) qeid.UseVHSPH(tru);
  if(!useVhsPh) qeid.UseVHSPH(fal);
  if(useQeEnuKinDif) qeid.UseQEkinEnuDif(tru);
  if(!useQeEnuKinDif) qeid.UseQEkinEnuDif(fal);
  if(useUserDefVariable) qeid.UseUserVar(tru);
  if(!useUserDefVariable) qeid.UseUserVar(fal);

  ///////////////////////////// main loop over snarls ////////////////////


  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;

    nevent = eventSummary->nevent;    
    run = ntpHeader->GetRun();
    snarl = ntpHeader->GetSnarl();


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

    detector = -1; 
    //get detector type:
    const DetectorType::Detector_t det 
      = ntpHeader->GetVldContext().GetDetector();
    if(det==DetectorType::kFar){
      detector = 2;
    }
    else if(det==DetectorType::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
    }

    if(nevent==0) continue;

    //fill PDFs with suitable events
    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;

      //////////////////////////////////////////////////////////////////////
      //zero all tree values

      flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0;

      mcevent = -1;
      is_fid = 0; is_cev = 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;

      reco_enu = 0; reco_emu = 0; reco_eshw = 0;  reco_wtd_eshw=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_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; 

      evtlength = trklength = shwlength=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;
      // kinematic quantities always positive
      // convention: -1 indicates unfilled variable (NC events)
      reco_y = reco_y = reco_Q2 = reco_W2 = -1;
     
      shw_nplanes_cut=shw_nplanes_raw=0;
      shw_ph_cut=shw_ph_raw=0.0;

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

      nduprs=0;
      dupphrs=0.;

      inu=0;

      ntotss=0;
      for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;}

      //////////////////////////////////////////////////////////////////////

      is_fid = InFidVol();

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

	inu = INu(mcevent);

      }

      if(PassBasicCuts(event)) {
	//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);

	pass_track = PassTrackCuts(track_index);
	//methods should all be protected against index = -1
	
	reco_emu          = RecoMKMuEnergy(emu_meth,track_index);

	int shower_index  = -1;

	// use the jim shower!
	LoadShower_Jim(i,shower_index,det);
	reco_eshw         = RecoShwEnergy(shower_index, det,1,!file_is_mc);
	if(shower_index>-1) reco_wtd_eshw = 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 is only useful for FD events
	  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));
	  
	  // angle reconstruction
	  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;
	    trkmomrange       = ntpTrack->momentum.range;
	    trkqp             = ntpTrack->momentum.qp;
	    if(trkqp!=0){
	      trkqp = 1.0/CorrectSignedMomentumFromCurvature(1.0/trkqp,
							     !file_is_mc, det);
	    }

	    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;
	    // 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
	    if(detector==1){//if it's the near detector
	      if(pitt_evt_class>0){ // if it was classified
		int do_meth=0; // emu reco method we should use
		if( (pitt_evt_class == 1) || (pitt_evt_class == 3) ) 
		  do_meth=2; // stoppers --> range
		else if( (pitt_evt_class == 2) || (pitt_evt_class == 4) ) 
		  do_meth=1; // punch-through --> curvature
		if(do_meth==1 || do_meth==2){
		   reco_emu = RecoMKMuEnergy(do_meth,track_index,
					     !file_is_mc,det);
		   reco_enu = reco_emu+reco_eshw;
		}
		emu_meth=do_meth; // was forgetting about this
	      }
	    }
	    else if(detector==2){//if it's the near detector
	      int do_meth=FarTrkContained(ntpTrack->end.x,
					  ntpTrack->end.y,
					  ntpTrack->end.z);
	      if(do_meth==1 || do_meth==2){
		reco_emu = RecoMKMuEnergy(do_meth,track_index);
		reco_enu = reco_emu+reco_eshw;
	      }
	      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);
	    }
	    
	    // pitt tracking quality cuts
	    if( (ntpTrack->fit.pass>0) 
		&& (ntpTrack->fit.chi2/ntpTrack->fit.ndof <20 )
		&& abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
	      pitt_trk_qual=1;
	    
	    ////////// parton model/DIS kinematics quantities///////// 
	    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);
	    ///////////////// END kinematics /////////////////////////

	    ///////// recalculate QE energy and Q2 variables /////////
	    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);
	    //////////////// END QE energy Q2 ////////////////////////

	  }

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

	  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;


	    ////////////////// access subshowers////////////
	    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
	    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();
	    
	    ///////// done with subshowers ////////////////
	  }

	  if(ntrack==1 && reco_enu>0.0 && reco_enu<140
	     && is_pitt_fid && pitt_trk_qual) pass=1;

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


      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);

      //runt event cut
      if((reco_enu-nvtp)/(reco_enu + nvtp)<0.95 && (nvtt*1e9)<50) continue;

      //track fit momentum error cut
      if(trkeqp!=0 && trkqp!=0){
	if(emu_meth==1 && abs(trkeqp/trkqp)>0.3) continue;
      }

      if(trkqp>0) continue;
      if(reco_trk_vtxz<0.5) continue;

      if(pass==1){
	qeid.FillPDFs(cc_nc,process,inu,this,i,reco_enu,reco_qe_enu,reco_emu,reco_W2,ntotss);
      }
    } // loop over events
  } // loop over entries
 
  //qeid.RebinPDFs(2);
  qeid.WritePDFs(QEtag);

}

//********************************************************
void MadMedQEAnalysis::CreatePDFCuts(Float_t flatEffVal, const char* QEtag)
{

  // is this MC??
  this->GetEntry(0);
  const bool file_is_mc
    =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);

  if(!file_is_mc) return;

  //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
  
  //Reco Quantities
  Int_t detector;         // Near = 1, Far = 2;
  Int_t run;              // run number
  Int_t snarl;            // snarl number
 
  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_cev;           // fully contained. 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)

  Int_t pass_track;       // the primary track passes track cuts
  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_eshw;      // reco shower energy  (shw.ph.gev/1.23)
  Float_t reco_mushw;     // showers along muon track  
  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_dircosneu; // reco track vtx dircos w.r.t neutrino
  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 shower vtx
  Float_t reco_shw_vtxy;      // reco y shower vtx
  Float_t reco_shw_vtxz;      // reco z shower 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 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

  // 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 end plane

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

  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;

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


  Float_t etu, etv, elu, elv;
  Float_t swu,swv,wswu,wswv; 

  // mark's QE id
  NewMadQEID qeid; 
  Bool_t fal = false;
  Bool_t tru = true;
  if(useRePhFracPdf) qeid.UseRePHfrac(tru);
  if(!useRePhFracPdf) qeid.UseRePHfrac(fal);
  if(useRecoW2Pdf) qeid.UseRecoWsqd(tru);
  if(!useRecoW2Pdf) qeid.UseRecoWsqd(fal);
  if(useNshowPdf) qeid.UseNshows(tru);
  if(!useNshowPdf) qeid.UseNshows(fal);
  if(useHpeakPdf) qeid.UseHoughPeak(tru);
  if(!useHpeakPdf) qeid.UseHoughPeak(fal);
  if(useNvhsHitsPdf) qeid.UseNVHShits(tru);
  if(!useNvhsHitsPdf) qeid.UseNVHShits(fal);
  if(useNSubS) qeid.UseNsubshows(tru);
  if(!useNSubS) qeid.UseNsubshows(fal);
  if(useVhsPh) qeid.UseVHSPH(tru);
  if(!useVhsPh) qeid.UseVHSPH(fal);
  if(useQeEnuKinDif) qeid.UseQEkinEnuDif(tru);
  if(!useQeEnuKinDif) qeid.UseQEkinEnuDif(fal);
  if(useUserDefVariable) qeid.UseUserVar(tru);
  if(!useUserDefVariable) qeid.UseUserVar(fal);

  qeid.ReadPDFs(QEtag);

  ///////////////////////////// main loop over snarls ////////////////////


  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;

    nevent = eventSummary->nevent;    
    run = ntpHeader->GetRun();
    snarl = ntpHeader->GetSnarl();


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

    detector = -1; 
    //get detector type:
    const DetectorType::Detector_t det 
      = ntpHeader->GetVldContext().GetDetector();
    if(det==DetectorType::kFar){
      detector = 2;
    }
    else if(det==DetectorType::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
    }

    if(nevent==0) continue;

    //fill PDFs with suitable events
    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;

      //////////////////////////////////////////////////////////////////////
      //zero all tree values

      flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0;

      mcevent = -1;
      is_fid = 0; is_cev = 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;

      reco_enu = 0; reco_emu = 0; reco_eshw = 0;  reco_wtd_eshw=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_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; 

      evtlength = trklength = shwlength=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;
      // kinematic quantities always positive
      // convention: -1 indicates unfilled variable (NC events)
      reco_y = reco_y = reco_Q2 = reco_W2 = -1;
      
      shw_nplanes_cut=shw_nplanes_raw=0;
      shw_ph_cut=shw_ph_raw=0.0;

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

      nduprs=0;
      dupphrs=0.;

      inu=0;

      ntotss=0;
      for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;}

      //////////////////////////////////////////////////////////////////////

      is_fid = InFidVol();

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

	inu = INu(mcevent);

      }

      if(PassBasicCuts(event)) {
	//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);

	pass_track = PassTrackCuts(track_index);
	//methods should all be protected against index = -1
	
	reco_emu          = RecoMKMuEnergy(emu_meth,track_index);

	int shower_index  = -1;

	// use the jim shower!
	LoadShower_Jim(i,shower_index,det);
	reco_eshw         = RecoShwEnergy(shower_index, det,1,!file_is_mc);
	if(shower_index>-1) reco_wtd_eshw = 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 is only useful for FD events
	  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));
	  
	  // angle reconstruction
	  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;
	    trkmomrange       = ntpTrack->momentum.range;
	    trkqp             = ntpTrack->momentum.qp;
	    if(trkqp!=0){
	      trkqp = 1.0/CorrectSignedMomentumFromCurvature(1.0/trkqp,
							     !file_is_mc, det);
	    }

	    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;
	    // 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
	    if(detector==1){//if it's the near detector
	      if(pitt_evt_class>0){ // if it was classified
		int do_meth=0; // emu reco method we should use
		if( (pitt_evt_class == 1) || (pitt_evt_class == 3) ) 
		  do_meth=2; // stoppers --> range
		else if( (pitt_evt_class == 2) || (pitt_evt_class == 4) ) 
		  do_meth=1; // punch-through --> curvature
		if(do_meth==1 || do_meth==2){
		   reco_emu = RecoMKMuEnergy(do_meth,track_index,
					     !file_is_mc,det);
		   reco_enu = reco_emu+reco_eshw;
		}
		emu_meth=do_meth; // was forgetting about this
	      }
	    }
	    else if(detector==2){//if it's the near detector
	      int do_meth=FarTrkContained(ntpTrack->end.x,
					  ntpTrack->end.y,
					  ntpTrack->end.z);
	      if(do_meth==1 || do_meth==2){
		reco_emu = RecoMKMuEnergy(do_meth,track_index);
		reco_enu = reco_emu+reco_eshw;
	      }
	      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);
	    }
	    
	    // pitt tracking quality cuts
	    if( (ntpTrack->fit.pass>0) 
		&& (ntpTrack->fit.chi2/ntpTrack->fit.ndof <20 )
		&& abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
	      pitt_trk_qual=1;
	    
	    ////////// parton model/DIS kinematics quantities///////// 
	    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);
	    ///////////////// END kinematics /////////////////////////

	    ///////// recalculate QE energy and Q2 variables /////////
	    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);
	    //////////////// END QE energy Q2 ////////////////////////

	  }

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

	  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;


	    ////////////////// access subshowers////////////
	    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
	    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();
	    
	    ///////// done with subshowers ////////////////
	  }

	  if(ntrack==1 && reco_enu>0.0 && reco_enu<140
	     && is_pitt_fid && pitt_trk_qual) pass=1;

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


      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);

      // DON'T PRE-SELECT EVENTS AS THIS AFFECTS THE FLATTENED EFFICIENCY IN THE FLUX EXTRACTION

      //runt event cut
      //if((reco_enu-nvtp)/(reco_enu + nvtp)<0.95 && (nvtt*1e9)<50) continue;

      //track fit momentum error cut
      //if(trkeqp!=0 && trkqp!=0){
      //  if(emu_meth==1 && abs(trkeqp/trkqp)>0.3) continue;
      //}

      //if(trkqp>0) continue;
      //if(reco_trk_vtxz<0.5) continue;

      //if(pass==1 && cc_nc==1 && process==1001){
      qeid.FillTrainingHists(this,i,reco_enu,reco_qe_enu,reco_emu,reco_W2,ntotss);
      //}
	
    } // loop over events
  } // loop over entries
 
  qeid.MakeCutHist(flatEffVal);
  qeid.WritePDFs(QEtag); 

}


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

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

//********************************************************
bool MadMedQEAnalysis::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()==DetectorType::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()==DetectorType::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);
    double radius = pow(vx-kXcenter,2) + pow(vy-kYcenter,2);

    radius=sqrt(radius);
    // restrict vertices to target region and 1m radius
    //if(vz>=kVetoEndZ && 
    //   vz<=kTargEndZ &&
    //   radius<=0.5) is_fid = true;    
    // 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 MadMedQEAnalysis::RecoMKMuEnergy(Int_t& opt,Int_t itrk, bool isdata, const DetectorType::Detector_t det){
  const float mumass=0.10566;
  if(LoadTrack(itrk)){
     float mr=ntpTrack->momentum.range;
     mr=CorrectMomentumFromRange(mr,isdata,det);
     float mc=0.0; // signed momentum from curvature
     if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
     mc=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 MadMedQEAnalysis::RecoShwEnergy(int opt, const DetectorType::Detector_t& det, int mode, bool isdata){
  
  Float_t result = 0;
  // Return Jim's calculation
  Bool_t ok = LoadShower(opt);
  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 = CorrectShowerEnergy(result,det,CandShowerHandle::kCC,mode,isdata);
  }
  return result;
  
}

//********************************************************
  /*
Float_t MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::PassBasicCuts(Int_t event){
  return kTRUE;
}
*/


//********************************************************
Float_t MadMedQEAnalysis::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 MadMKAnalysys::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 MadMedQEAnalysis::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(uint 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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::SetBECFile(const char* f){
  if(f){
    fBECConfig.Set("beam:flux_file",f);
  }
}


//********************************************************
Int_t MadMedQEAnalysis::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<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 MadMedQEAnalysis::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];
      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 MadMedQEAnalysis::GetEvtTimeChar(double &evttimemin, double &evttimemax, double &evttimeln)
{
  //cut and pasted, with minor changes from MadDpAnalysis.cxx
  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()==DetectorType::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      LoadStrip(stp_index); 
      Double_t striptime=ntpStrip->time1-trgtime;
      if(striptime<=timemin) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
  }	 
  
  if(ntpHeader->GetVldContext().GetDetector()==DetectorType::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      LoadStrip(stp_index);
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      Double_t striptime;
      if(striptime1>0 && striptime0<0)  striptime=striptime1;
      if(striptime0>0 && striptime1<0)  striptime=striptime0;
      if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
      if(striptime<=timemin) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
  }

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

}

//********************************************************
void MadMedQEAnalysis::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()==DetectorType::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      LoadStrip(stp_index); 
      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()==DetectorType::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      LoadStrip(stp_index);
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      Double_t striptime;
      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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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 MadMedQEAnalysis::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];
    LoadStrip(si);
    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);
}