MadMKAnalysis.cxx

Go to the documentation of this file.

//
// MadMKAnalysis
//
// Mike's MadAnalysis.  Mostly I just make a simple ntuple with it.
//
// Created:  M. Kordosky -- December, 2004
//
// $Author: rhatcher $ 
//
// $Revision: 1.61 $
// 
// $Name:  $
//
// $Id: MadMKAnalysis.cxx,v 1.61 2008/11/19 17:16:21 rhatcher Exp $
//
// $Log: MadMKAnalysis.cxx,v $
// Revision 1.61  2008/11/19 17:16:21  rhatcher
// replace non-standard "uint" with "unsigned int".
//
// Revision 1.60  2007/08/23 15:58:45  kordosky
// small updates, not committed before
//
// Revision 1.59  2007/06/10 19:29:12  kordosky
// reweight to change uniform density to WS
//
// Revision 1.58  2007/06/05 17:10:41  kordosky
// blah
//
// Revision 1.57  2007/05/08 12:50:48  kordosky
// support intranuke reweighting
//
// Revision 1.56  2007/02/01 22:01:19  rhatcher
// Changes to ROOT headers means we need an explicit #include "TClonesArray.h"
//
// Revision 1.55  2007/02/01 17:48:17  kordosky
// New version of EnergyCorrections. EnergyCorrections is now a namespace and calling code should invoke "using namespace EnergyCorrections". This is not expected
// to cause problems since in the previous version all the correction functions were at global scope. Indeed, code including EnergyCorrections.h have been modified. Will handle cedar and birch by using static variables fVersion and fSubVersion, modified via SetCorrectionVersion(). Users need to do this, perhaps in their job macro. So far there is only one correction of each type and fSubVersion should be held at the default 0. Also, note the utility VersionFromFilename().
//
// Revision 1.54  2007/01/30 03:17:21  rhatcher
// No longer user BeamType "class" found in MCReweight, but rather the
// BeamType "namespace" found in Conventions package.
// Also covert to using "Detector::" rather than "DetectorType::", though
// for now they are synonyms.
//
// Revision 1.53  2006/10/16 19:53:30  kordosky
// cope with CEDAR strip indices set to -1. allow user to force beam type (useful for 150cm data). Open files with TFile::Open() to allow natural dcache access. IMPORTANT: Big change to pans. Tacking on a spill marker, looks like event=-1 with no muon or shower energy, pass and is_fid set to fail. Beam variables and nevent are correct for these entries. Modification intended to ease POT counting
//
// Revision 1.52  2006/05/22 19:16:38  rhatcher
// Initialize some local variables that gcc worries might be used before
// they are set.
//
// Revision 1.51  2006/03/10 18:57:21  kordosky
// correct MadDpID and MadMKAnalysis to account for 1.8% shift in pulseheight related PID quantities
//
// Revision 1.50  2006/03/09 19:33:30  kordosky
// Properly do energy corrections in MK and DpAnalysis
//
// Revision 1.49  2006/03/08 19:32:54  kordosky
// Account for detector and data/mc dependent energy corrections in MadQuantities. Attempt to fix DpAnalysis to work with the new changes.
//
// Revision 1.48  2006/02/16 19:26:09  kordosky
// upgrades
//
// Revision 1.47  2006/02/01 20:17:54  kordosky
// add shower truth/ completeness.
//
// Revision 1.46  2006/01/27 17:04:17  vahle
// Adding flux information to MadMKAnalysis code
//
// Revision 1.45  2006/01/27 11:18:22  kordosky
// upgrades
//
// Revision 1.44  2006/01/26 23:41:34  kordosky
// upgrades.
//
// Revision 1.43  2006/01/20 17:30:04  kordosky
// Enable user control of shower marker color and style.
//
// Revision 1.42  2005/12/15 21:56:36  kordosky
// fix to use datautil version of energycorrections.h
//
// Revision 1.41  2005/12/15 02:42:43  kordosky
// tweaks to class constructors
//
// Revision 1.40  2005/12/14 20:35:59  kordosky
// Updated energy corrections.
//
// Revision 1.39  2005/11/18 17:13:16  vahle
// Add RemoteSpillType to MadMKAnalysis
//
// Revision 1.38  2005/11/18 16:26:45  vahle
// Add RemoteSpillType to MadMKAnalysis
//
// Revision 1.37  2005/11/17 16:59:35  vahle
// Must set time diff AFTER setting spill
//
// Revision 1.36  2005/11/17 16:34:14  vahle
// Must set time diff manually when trying to use BMSpillAna with data base instead of ntuple beam records
//
// Revision 1.35  2005/11/17 15:50:22  vahle
// Added checks for track vertex in fid volume.  Switched goodbeam varible to use BMSpillAna
//
// Revision 1.34  2005/11/17 13:21:58  vahle
// Added a count of the number of reconstructed strips and the ph of reconstructed strips that show up in another event
//
// Revision 1.33  2005/11/14 12:05:00  vahle
// Adding Brians lowph ratio variable for cutting fd junk events to MadMKAnalysis ntuple
//
// Revision 1.32  2005/11/04 15:27:45  vahle
// Remove all vestiges of old beam monitoring ntuple lookup, change the way track containment is decided in the FD
//
// Revision 1.31  2005/10/19 23:56:42  kordosky
// *** empty log message ***
//
// Revision 1.30  2005/10/19 20:36:40  vahle
// Fix pot counting in FD
//
// Revision 1.29  2005/10/19 19:55:37  vahle
// Added new variables to MadMKAnalysis ntuple, namely li time, demux status; fixed various small bugs
//
// Revision 1.28  2005/10/07 18:33:00  vahle
// Added late activity variables to MadMKAnalysis
//
// Revision 1.27  2005/10/07 18:05:28  kordosky
// a few changes to analysis, refinement of pdfs. MadBase not loading mcHeader correctly when GetEntry() called.  Suprisingly, it had no effect up until now.
//
// Revision 1.26  2005/10/06 18:18:09  vahle
// bug fix, use plane number instead of z when asking if strip is within a plane window
//
// Revision 1.25  2005/10/06 17:52:08  vahle
// Added early activity variables to ntuple
//
// Revision 1.24  2005/10/06 17:29:28  vahle
// Added early activity variables to MadMKAnalysis
//
// Revision 1.23  2005/10/06 15:03:26  kordosky
// various improvements, pids finally integrated into pan construction routine, still not fully debugged though
//
// Revision 1.22  2005/09/23 11:48:22  kordosky
//  a few additions for R1.18. code meant to handle R1.16 and R1.18 ntpshowers not working properly. only handles R1.18 correctly.
//
// Revision 1.21  2005/09/14 13:32:45  kordosky
// DpID and NsID classes now fully working and validated with MadTestAnalysis by comparing nannpid vs. annpid and npid0 vs. pid0.  pdf files and weights stored in data subdirectory. FD pdfs needed for DpID.
//
//

#include <set>
#include <vector>
#include "TClonesArray.h"
#include "TSystem.h"
#include "TRandom.h"
#include "TVectorD.h"

#include "Mad/MadMKAnalysis.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 "NeugenInterface/neugen_config.h"
#include "NeugenInterface/neugen_wrapper.h"
#include "NeugenInterface/inuke_rw.h"

#include "Conventions/BeamType.h"

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

using namespace EnergyCorrections;

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

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

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

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

}

//********************************************************
MadMKAnalysis::MadMKAnalysis(JobC *j,string path,int entries) : MadAnalysis(j, path, entries)
{
  fMCBeam=BeamType::kUnknown;
  fDataUseMCBeam=false;
  fDoInukeReweight=false;
  fNRandomSets=0;
  fOnlyWriteFidEvents=false;
  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  Clear();
  isMC = true;
  isTH = true;
  isEM = true;
  Nentries = entries;
  whichSource = 1;
  jcPath = path;
  fJC = j;
  fChain = NULL;


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

}

//********************************************************
MadMKAnalysis::~MadMKAnalysis()
{

}

void MadMKAnalysis::CreatePAN(std::string tag, const char* /*bmonpath*/)
{
  
  // 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: MadMKAnalysis::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_mushw;     // showers along muon track
  Float_t reco_eshw;      // reco shower energy (largest shower
  Float_t reco_wtd_eshw;    // as above but weighted 
  Float_t reco_vtx_eshw;    // reco shower energy (=0 if no vertex shower)
  Float_t reco_vtx_wtd_eshw;// as above but weighted (=0 if no vertex shower)
  
  Float_t reco_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 track vtx
  Float_t reco_shw_vtxy;      // reco y track vtx
  Float_t reco_shw_vtxz;      // reco z track vtx
  Float_t reco_shw_vtxr;// radial position of vertex


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

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


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

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

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


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

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



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

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

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

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


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

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

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

  //time last li event in snarl, -1 if none;
  double litime;

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

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


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

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

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

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

  //Reco
  tree->Branch("detector",&detector,"detector/I",32000);
  tree->Branch("run",&run,"run/I",32000);
  tree->Branch("snarl",&snarl,"snarl/I",32000);
  tree->Branch("subrun",&subrun,"subrun/I",32000);
  tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000);
  tree->Branch("trgtime",&trgtime,"trgtime/D",32000);
  tree->Branch("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 med_qe_pid;
  int med_qe_pass;
  tree->Branch("med_qe_pid",&med_qe_pid,"med_qe_pid/F",32000);
  tree->Branch("med_qe_pass",&med_qe_pass,"med_qe_pass/I",32000);

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


  tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
  tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
  tree->Branch("reco_mushw",&reco_mushw,"reco_mushw/F",32000);
  tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
  tree->Branch("reco_wtd_eshw",&reco_wtd_eshw,"reco_wtd_eshw/F",32000);
  tree->Branch("reco_vtx_eshw",&reco_vtx_eshw,"reco_vtx_eshw/F",32000);
  tree->Branch("reco_vtx_wtd_eshw",&reco_vtx_wtd_eshw,"reco_vtx_wtd_eshw/F",32000);
  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("shw_nplanes_cut",&shw_nplanes_cut,"shw_nplanes_cut/I",32000);
  tree->Branch("shw_nplanes_raw",&shw_nplanes_raw,"shw_nplanes_raw/I",32000);

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


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

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

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

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


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


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

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

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

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

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

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

  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 a unsigned 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;


  TTree* inukerw = 0;
  int n_inuke_wts=0;
  float* inuke_wts=0;
  double nucrad=-1.0;
  double wrad=1.0;
  std::vector<inuke_reweight::parameter_set> psets;
  
  if(file_is_mc && fDoInukeReweight){
    
    inuke_reweight::parameter_set sigmas;

    neugen_config nc("MODBYRS");
    nc.set_config_no(4);
    neugen_wrapper nw(&nc);
    nw.print_configuration();
    gSystem->Setenv("NEUGEN3_XSEC","MODBYRS4_xsec.dat");
    nw.begin_generation();

    /*    
    sigmas.pi_fate.inel=0.2;
    sigmas.pi_fate.abs=0.25;
    sigmas.pi_scale.ft=0.5;
    sigmas.pi_scale.xsec=0.1;    
    //    sigmas.pn_scale.ft=0.5;
    sigmas.pn_scale.xsec=0.1;
    */

    sigmas.pi_fate.elas=0.1;
    sigmas.pi_fate.inel=0.4;
    sigmas.pi_fate.abs=0.3;
    sigmas.pi_fate.cex=0.49;
    sigmas.pi_fate.piprod=0.2;

    sigmas.pi_scale.xsec=0.1;

    sigmas.pn_fate.abs=0.2;
    sigmas.pn_fate.piprod=0.2;

    sigmas.pn_scale.xsec=0.1;

    sigmas.pi_scale.ft=0.5; // apply to both pi and p

    std::cout<<"Inuke Reweight sigmas: "<<std::endl;
    sigmas.print(std::cout);

    if(fNRandomSets<=0) {
      inuke_reweight::generate_1sigma_shifts(sigmas,psets);    
      for (unsigned int iset=0; iset<psets.size(); iset++){
        psets[iset].pi_fate.set_abs_fates();
        psets[iset].pn_fate.set_abs_fates();
        psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime
        // correlate change in pn xsec with change in pn elastic fate
        if(psets[iset].pn_scale.xsec!=0) {
          psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec;
        }
      }
    }
    else{
      inuke_reweight::parameter_limits limits(sigmas,2.0);
      limits.upper.pi_scale.ft=0.5;  // special case, 1 sigma, but uniform dist
      limits.lower.pi_scale.ft=-0.5;
      std::cout<<"Inuke Reweight Limits: "<<std::endl;
      limits.print(std::cout);
      gRandom->SetSeed(1);
      for(int i=0; i<fNRandomSets; i++){
        inuke_reweight::parameter_set p;
        inuke_reweight::generate_uncor_shifts(sigmas,limits,p);
        psets.push_back(p);
      }
      // massage parameters
      for(int iset=0; iset<fNRandomSets; iset++){
        psets[iset].pi_fate.set_abs_fates();
        psets[iset].pn_fate.set_abs_fates();
        psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime
        // correlate change in pn xsec with change in pn elastic fate
        if(psets[iset].pn_scale.xsec!=0) {
          psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec;
        }
      }


    }

    n_inuke_wts=psets.size();
    inuke_wts = new float[n_inuke_wts];
    
    //    inukerw = new TTree("inukerw","intranuke weights");
    inukerw=tree;
    inukerw->Branch("r",&run,"r/I");
    inukerw->Branch("s",&snarl,"s/I");
    inukerw->Branch("e",&event,"e/I");
    inukerw->Branch("nwts",&n_inuke_wts,"nwts/I");
    inukerw->Branch("wts",inuke_wts,"wts[nwts]/F");          
    inukerw->Branch("nucrad",&nucrad,"nucrad/D");
    inukerw->Branch("wrad",&wrad,"wrad/D");

    // save psets
    // 9*2 fates + 3*2 scales = 24
    for(int ipar=0; ipar<24; ipar++){
      TVectorD* vptr = new TVectorD(n_inuke_wts);
      inukerw->GetUserInfo()->Add(vptr);
    }
    for(int iset=0; iset<n_inuke_wts; iset++){
      for(int ipar=0; ipar<24; ipar++){ 
        TVectorD* vptr = 
          dynamic_cast<TVectorD*>(inukerw->GetUserInfo()->At(ipar));    
        (*vptr)[iset]=psets[iset].get_par(ipar);
      }      
    }


  }



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

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

  // mark's QE id
  //  MadQEID qeid;
  //  qeid.ReadPDFs("/afs/fnal.gov/files/data/minos/d04/kordosky/madqe_files/QEpdfs.root");

  
  //  bool beamdb=false;
  //  bool checkeddb=false;
  int lastrun=-1;
  for(int ii=0;ii<Nentries;ii++){   
    if(ii%1000==0) std::cout << float(ii)*100./float(Nentries) 
                            << "% done" << std::endl;

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

    /*
    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 Detector::Detector_t det 
      = ntpHeader->GetVldContext().GetDetector();
    if(det==Detector::kFar){
      detector = 2;
    }
    else if(det==Detector::kNear){
      detector = 1;     
    }
    // deal with beam type as best we can
    BeamType::BeamType_t current_beam=BeamType::kUnknown;
    if(file_is_mc){
      // for MC set beam to fMCbeam
      current_beam=fMCBeam;
      // for data we will try to deduce from beam monitoring
    }


    // first thing, if data do beam monitoring
    hbw=vbw=hpos1=vpos1=hpos2=vpos2=hornI=closestspill=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;
      }


      if(fDataUseMCBeam==true){
        //allow us to force the data to use the MC beam
        current_beam=fMCBeam;
      }
      else{
        // attempt to deduce beam type from beam monitoring
        current_beam = BeamType::FromBeamMon(beamcnf);
      }
      
      if(goodbeam==1&&spilltype!=3){//don't count fake spills
        pot+=tortgt;
        pottor101+=tor101;
        pottortgt+=tortgt;
      }
    }

    // Write a special record for each spill!
    event=-1; ntrack=-1; nshower=-1; reco_emu=0.0; reco_enu=0.0; reco_eshw=0;
    true_enu=0.0; true_emu=0.0; true_eshw=0.0;
    pass=0; is_fid=0;
    tree->Fill();
    for(int iwt=0; iwt<n_inuke_wts; iwt++) inuke_wts[iwt]=-1.0;
    if(inukerw!=0 && inukerw!=tree) inukerw->Fill();

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

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

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

      fluxrun=0; fluxevtno=0; nudxdz=0.; nudydz=0.; nupz=0.; nuenergy=0.;
      nudxdznear=0.; nudydznear=0.; nuenergynear=0.; nwtnear=0.; nudxdzfar=0.;
      nudydzfar=0.; nuenergyfar=0.; nwtfar=0.; norig=0; ndecay=0; ntype=0;
      vx=0.; vy=0.; vz=0.; pdpx=0.; pdpy=0.; pdpz=0.; ppdxdz=0.; ppdydz=0.; pppz=0.;
      ppenergy=0.; ppmedium=0; ptype=0; ppvx=0.; ppvy=0.; ppvz=0.;
      muparpx=0.; muparpy=0.; muparpz=0.; mupare=0.; necm=0.; nimpwt=0.;
      xpoint=0.; ypoint=0.; zpoint=0.; tvx=0.; tvy=0.; tvz=0.; tpx=0.;
      tpy=0.; tpz=0.; tptype=0; tgen=0;
      
      flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0;

      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_x = 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;}

      med_qe_pid=-1001.0;
      med_qe_pass=0;
      niki_cc_pid=-999;
      dave_cc_pid=-999;

      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;
            trkmomrange = CorrectMomentumFromRange(trkmomrange,!file_is_mc,det);
            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);
            
            // CEDAR apparently has 5 plane long tracks, with ndf==0
            Double_t tmp_ndof=ntpTrack->fit.ndof;
            if(tmp_ndof<=0.0) tmp_ndof=1.0;
            // pitt tracking quality cuts
            if( (ntpTrack->fit.pass>0) 
                && (ntpTrack->fit.chi2/tmp_ndof <20 )
                && abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
              pitt_trk_qual=1;
            
            const double M=(0.93827 + 0.93957)/2.0; // <nucleon mass>
            if(reco_emu>0) reco_Q2 = 
                             2*reco_enu*reco_emu*(1.0 - reco_dircosneu);
            reco_W2 = M*M - reco_Q2 + 2*M*reco_eshw;      
            if(reco_enu>0) reco_y = reco_eshw/reco_enu;
            if(reco_eshw>0 && reco_Q2>0) reco_x =  reco_Q2/(2*M*reco_eshw);

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

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

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

            if(det==Detector::kFar && !file_is_mc){
              dpid.SetPHCorrection(1.018);
            }
            if(dpid.ChoosePDFs(det,current_beam)){
              dave_cc_pid = dpid.CalcPID(this,event,0);
            }

            TObject *recobj=0;
            if(record!=0){
              recobj=record;
            }
            else{
              recobj=strecord;
            }
              
    
            LoadTrack(track_index);
            //figure out planes in trk > 1.5 pe
            float temp_ph;
            trk_nplanes=NPlanesInObj(recobj,ntpTrack,1.5,temp_ph);
            //figure out strips in trk > 1.5 pe
            trk_nstrips=NStripsInObj(recobj,ntpTrack,1.5,temp_ph);

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

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

          float temp_ph;
          //figure out planes in slice > 1.5 pe
          slc_nplanes=NPlanesInObj(recobj,ntpSlice,1.5,temp_ph);
          //figure out planes in evt > 1.5 pe
          evt_nplanes=NPlanesInObj(recobj,ntpEvent,1.5,temp_ph);
          //figure out strips in slice > 1.5 pe
          slc_nstrips=NStripsInObj(recobj,ntpSlice,1.5,temp_ph);
          //figure out strips in evt > 1.5 pe
          evt_nstrips=NStripsInObj(recobj,ntpEvent,1.5,temp_ph);
          
          
          if(shower_index!=-1 && reco_eshw>0){ // case of no vertex shower
            LoadShower(shower_index);
            shwlength=ntpShower->plane.n;
            shw_nplanes_raw=NPlanesInObj(recobj,ntpShower,0.0,shw_ph_raw);
            shw_nplanes_cut=NPlanesInObj(recobj,ntpShower,1.5,shw_ph_cut);
            
            shwend=ntpShower->plane.end;
            shwbeg=ntpShower->plane.beg;
            reco_shw_vtxx=ntpShower->vtx.x;
            reco_shw_vtxy=ntpShower->vtx.y;
            reco_shw_vtxz=ntpShower->vtx.z;
            
            
            int ncl = ntpShower->ncluster;

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

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

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

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

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

      if(file_is_mc){
        if(gnumi->Status()) {
          if(isST) gnumi->GetParent(strecord,mcevent,*nuparent);
          else gnumi->GetParent(mcrecord,mcevent,*nuparent);
        }               
      }

      if(fDoInukeReweight && file_is_mc && inukerw!=0) {
        
        if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents){
          std::vector<double> weights(n_inuke_wts,1.0);
          //    std::cout<<"snarl, event: "<<snarl<<", "<<event<<std::endl;
          inuke_reweight::calc_weights(strecord,ntpTruth,psets,weights,nucrad,wrad,false);
          for (unsigned int iwt=0; iwt<weights.size(); iwt++) inuke_wts[iwt]=weights[iwt];
        }
        // nuclear radius and weighting factor

      }

      if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents) {
        tree->Fill(); if(inukerw!=0 && inukerw!=tree) inukerw->Fill();
      }


    } // loop over events
  } // loop over entries
  
  // delete nuparent; // should I do this?
  save.cd();
  pottree->Fill();
  pottree->Write();
  tree->Write();if(inukerw!=0 && inukerw!=tree) inukerw->Write();

  save.Write();
  save.Close();

}

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

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

bool MadMKAnalysis::InFidVol(float vx, float vy, float vz){
  // assumes event has been loaded
  // check that vertex is in fiducial volume
  bool is_fid=false;

  // FD: from doc-1351-v2 : there is no coil cut!
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
        //fiducial criteria on vtx for far detector
    is_fid = true;
    if(vz<0.5 || vz>28.0 ||   //ends
       (vz<16.2&&vz>14.3) ||  //between SMs
       ((vx*vx)+(vy*vy))>14.0) is_fid = false;
  }
  else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
    //fiducial criteria on vtx for near detector

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

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


    is_fid = false;
    // MAK: Sept 19, 2005: Modified to account for beam slope
    //    double radius = pow(vx-kXcenter,2) 
    //      + pow(vy-(kYcenter -vz*TMath::Tan(0.058)),2);
    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 MadMKAnalysis::RecoMKMuEnergy(Int_t& opt,Int_t itrk, bool isdata, const Detector::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 MadMKAnalysis::RecoShwEnergy(int opt, const Detector::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::PassBasicCuts(Int_t event){
  return kTRUE;

}
*/

Float_t MadMKAnalysis::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 MadMKAnalysis::NDRadialFidVol(Float_t x,Float_t y, Float_t z){

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

Bool_t MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::SetBECFile(const char* f){
  if(f){
    fBECConfig.Set("beam:flux_file",f);
  }
}

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

  std::set<UShort_t> planes;

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

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

  if(slc!=0){
    stripindex=slc->stp;
    nstrip=slc->nstrip;
    //    cout<<"Counting planes in slice"<<endl;
  }
  else if(trk!=0){
    stripindex=trk->stp;
    nstrip=trk->nstrip;
    //    cout<<"Counting planes in trk"<<endl;
  }
  else if(evt!=0){
    stripindex=evt->stp;
    nstrip=evt->nstrip;
    //    cout<<"Counting planes in event"<<endl;
  }
  else if(shw!=0){
    stripindex=shw->stp;
    nstrip=shw->nstrip;
    //    cout<<"Counting planes in event"<<endl;
  }
  
  TClonesArray *striplist=0; 
  int TOTSTRIPS=0;
  if(srrec!=0){
    striplist = srrec->stp;
    TOTSTRIPS=srrec->evthdr.nstrip;
    //    cout<<"Found a sr"<<endl;
  }
  else{
    striplist = strec->stp;
    TOTSTRIPS=strec->evthdr.nstrip;
    //    cout<<"found a st"<<endl;
  }
  
  if(striplist!=0){
    //    cout<<"Looping over strip list, nstrip="<<nstrip<<endl;
    for(int i=0;i<nstrip;i++){
      int sindex = stripindex[i];
      if(sindex<0) continue;
      if(sindex<TOTSTRIPS){       
        NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>
          ((*striplist)[sindex]);
        if(!strip) continue;
        //      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]);
      if(!strip) continue;
      //      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 MadMKAnalysis::NStripsInObj(TObject *rec, TObject *obj, float phcut, float& sumph)
{
  sumph=0.0;
  std::vector<UShort_t> strips;

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

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

  if(slc!=0){
    stripindex=slc->stp;
    nstrip=slc->nstrip;
  }
  else if(trk!=0){
    stripindex=trk->stp;
    nstrip=trk->nstrip;
  }
  else if(evt!=0){
    stripindex=evt->stp;
    nstrip=evt->nstrip;
  }
  else if(shw!=0){
    stripindex=shw->stp;
    nstrip=shw->nstrip;
  }
  
  TClonesArray *striplist=0; 
  int TOTSTRIPS;
  if(srrec!=0){
    striplist = srrec->stp;
    TOTSTRIPS=srrec->evthdr.nstrip;
  }
  else{
    striplist = strec->stp;
    TOTSTRIPS=strec->evthdr.nstrip;
  }
  
  //  cout<<"CHECKING nstrips="<<nstrip<<endl;
  if(striplist!=0){
    for(int i=0;i<nstrip;i++){
      int sindex=stripindex[i];
      if(sindex<0) continue;
      NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[sindex]);
      if(!strip) continue;
      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) continue;
      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 MadMKAnalysis::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()==Detector::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime=ntpStrip->time1-trgtime;
      if(striptime<=timemin) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
  }      
  
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      Double_t striptime = 0;
      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 MadMKAnalysis::GetEvtTimeCharPHC(double &evttimemin, double &evttimemax, double &evttimeln)
{
  //same as GetEvtTimeChar but only looking at hits with 
  //pulse height > 200 sig cor (~2pe) in near det
  //pulse height > 160 sig cor in far det
  double trgtime=eventSummary->trigtime;
  double timemax=0.;
  double timemin=1.e10;

  //  Find max min time of events by loading strips for ND 
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime=ntpStrip->time1-trgtime;
      if(ntpStrip->ph1.sigcor>200){
        if(striptime<=timemin) timemin=striptime;
        if(striptime>=timemax) timemax=striptime;
      }
    }
  }      
  
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
    for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
      Int_t stp_index=((ntpEvent->stp)[evss]);
      if(!LoadStrip(stp_index)) continue;
      Double_t striptime1=ntpStrip->time1-trgtime;
      Double_t striptime0=ntpStrip->time0-trgtime;
      Double_t striptime = 0;
      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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::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 MadMKAnalysis::DupRecoStrips(int evidx, int &nduprs, float &dupphrs)
{
  //zero the variables
  nduprs=0;
  dupphrs=0.;

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

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

  //reload the current event
  LoadEvent(evidx);
}

---