#define MedPanAna_cxx
#include "MedPanAna.h"
#include <TH1.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <iostream>
#include <cstdlib>
#include <cmath>

using namespace std;

//***************************************************
void MedPanAna::AnalyseMC(Bool_t selectCC, Int_t processToKeep, Bool_t lookAtNC, const char* tag)
{

  if(fChain==0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  Int_t nBytes = 0, nb = 0;

  TFile* fout = new TFile(tag,"recreate");
  TDirectory* dir = gDirectory;

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

  // Declare histograms

  // Histogram array convention:
  //   [0] - all events
  //   [1] - 'safe' Dis    Q2>1,  W2>4
  //   [2] - 'low q2' Dis  Q2<=1, W2>4
  //   [3] - QE/Delta      W2<1.69
  //   [4] - 'transition'  W2>=1.69 && W2<=4

  // 1D

  TH1F* enuB4[5];
  TH1F* enuB4Rw[5];
  TH1F* enuAft[5];
  TH1F* enuAftRw[5];

  TH1F* emuB4[5];
  TH1F* emuB4Rw[5];
  TH1F* emuAft[5];
  TH1F* emuAftRw[5];

  TH1F* eshwB4[5];
  TH1F* eshwB4Rw[5];
  TH1F* eshwAft[5];
  TH1F* eshwAftRw[5];

  TH1F* enuqeB4[5];
  TH1F* enuqeB4Rw[5];
  TH1F* enuqeAft[5];
  TH1F* enuqeAftRw[5];

  TH1F* thetamuB4[5];
  TH1F* thetamuB4Rw[5];
  TH1F* thetamuAft[5];
  TH1F* thetamuAftRw[5];

  TH1F* q2qeB4[5];
  TH1F* q2qeB4Rw[5];
  TH1F* q2qeAft[5];
  TH1F* q2qeAftRw[5];

  TH1F* q2B4[5];
  TH1F* q2B4Rw[5];
  TH1F* q2Aft[5];
  TH1F* q2AftRw[5];

  TH1F* xB4[5];
  TH1F* xB4Rw[5];
  TH1F* xAft[5];
  TH1F* xAftRw[5];

  TH1F* yB4[5];
  TH1F* yB4Rw[5];
  TH1F* yAft[5];
  TH1F* yAftRw[5];

  TH1F* w2B4[5];
  TH1F* w2B4Rw[5];
  TH1F* w2Aft[5];
  TH1F* w2AftRw[5];

  // 2D

  TH2F* thetamuVemuB4[5];
  TH2F* thetamuVemuB4Rw[5];
  TH2F* thetamuVemuAft[5];
  TH2F* thetamuVemuAftRw[5];

  TH2F* thetamuVeshwB4[5];
  TH2F* thetamuVeshwB4Rw[5];
  TH2F* thetamuVeshwAft[5];
  TH2F* thetamuVeshwAftRw[5];

  TH2F* yVenuB4[5];
  TH2F* yVenuB4Rw[5];
  TH2F* yVenuAft[5];
  TH2F* yVenuAftRw[5];

  TH2F* q2VenuB4[5];
  TH2F* q2VenuB4Rw[5];
  TH2F* q2VenuAft[5];
  TH2F* q2VenuAftRw[5];

  TH2F* w2VenuB4[5];
  TH2F* w2VenuB4Rw[5];
  TH2F* w2VenuAft[5];
  TH2F* w2VenuAftRw[5];

  TH2F* xVenuB4[5];
  TH2F* xVenuB4Rw[5];
  TH2F* xVenuAft[5];
  TH2F* xVenuAftRw[5];

  TH2F* xVq2B4[5];
  TH2F* xVq2B4Rw[5];
  TH2F* xVq2Aft[5];
  TH2F* xVq2AftRw[5];

  TH2F* rmtotenuVenuB4[5];
  TH2F* rmtotenuVenuB4Rw[5];
  TH2F* rmtotenuVenuAft[5];
  TH2F* rmtotenuVenuAftRw[5];

  TH2F* rmtoteshwVeshwB4[5];
  TH2F* rmtoteshwVeshwB4Rw[5];
  TH2F* rmtoteshwVeshwAft[5];
  TH2F* rmtoteshwVeshwAftRw[5];

  TH2F* rmtotemuVemuB4[5];
  TH2F* rmtotemuVemuB4Rw[5];
  TH2F* rmtotemuVemuAft[5];
  TH2F* rmtotemuVemuAftRw[5];

  TH2F* rmtotthetamuVthetamuB4[5];
  TH2F* rmtotthetamuVthetamuB4Rw[5];
  TH2F* rmtotthetamuVthetamuAft[5];
  TH2F* rmtotthetamuVthetamuAftRw[5];

  TH2F* rmtotthetamuVemuB4[5];
  TH2F* rmtotthetamuVemuB4Rw[5];
  TH2F* rmtotthetamuVemuAft[5];
  TH2F* rmtotthetamuVemuAftRw[5];

  TH2F* enuVeshwB4[5];
  TH2F* enuVeshwB4Rw[5];
  TH2F* enuVeshwAft[5];
  TH2F* enuVeshwAftRw[5];

  TH2F* emuVeshwB4[5];
  TH2F* emuVeshwB4Rw[5];
  TH2F* emuVeshwAft[5];
  TH2F* emuVeshwAftRw[5];

  TH2F* enuqeVeshwB4[5];
  TH2F* enuqeVeshwB4Rw[5];
  TH2F* enuqeVeshwAft[5];
  TH2F* enuqeVeshwAftRw[5];

  TH2F* q2qeVeshwB4[5];
  TH2F* q2qeVeshwB4Rw[5];
  TH2F* q2qeVeshwAft[5];
  TH2F* q2qeVeshwAftRw[5];

  TH2F* q2VeshwB4[5];
  TH2F* q2VeshwB4Rw[5];
  TH2F* q2VeshwAft[5];
  TH2F* q2VeshwAftRw[5];

  TH2F* xVeshwB4[5];
  TH2F* xVeshwB4Rw[5];
  TH2F* xVeshwAft[5];
  TH2F* xVeshwAftRw[5];

  TH2F* yVeshwB4[5];
  TH2F* yVeshwB4Rw[5];
  TH2F* yVeshwAft[5];
  TH2F* yVeshwAftRw[5];

  TH2F* w2VeshwB4[5];
  TH2F* w2VeshwB4Rw[5];
  TH2F* w2VeshwAft[5];
  TH2F* w2VeshwAftRw[5];

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

  TH2F* trkVtxB4 = new TH2F("trkVtxB4","trkVtxB4",800,0.5,2.5,600,-1.5,1.5);
  TH2F* trkVtxAft = new TH2F("trkVtxAft","trkVtxAft",800,0.5,2.5,600,-1.5,1.5);

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

  // Initialise histograms

  char name[50];
  for(Int_t i=0;i<5;i++){
    
    // 1D
    
    sprintf(name,"enuB4[%d]",i);
    enuB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuB4Rw[%d]",i);
    enuB4Rw[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuAft[%d]",i);
    enuAft[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuAftRw[%d]",i);
    enuAftRw[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"emuB4[%d]",i);
    emuB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"emuB4Rw[%d]",i);
    emuB4Rw[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"emuAft[%d]",i);
    emuAft[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"emuAftRw[%d]",i);
    emuAftRw[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"eshwB4[%d]",i);
    eshwB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"eshwB4Rw[%d]",i);
    eshwB4Rw[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"eshwAft[%d]",i);
    eshwAft[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"eshwAftRw[%d]",i);
    eshwAftRw[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"enuqeB4[%d]",i);
    enuqeB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuqeB4Rw[%d]",i);
    enuqeB4Rw[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuqeAft[%d]",i);
    enuqeAft[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuqeAftRw[%d]",i);
    enuqeAftRw[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"thetamuB4[%d]",i);
    thetamuB4[i] = new TH1F(name,name,200,-0.1,1.9);
    sprintf(name,"thetamuB4Rw[%d]",i);
    thetamuB4Rw[i] = new TH1F(name,name,200,-0.1,1.9);
    sprintf(name,"thetamuAft[%d]",i);
    thetamuAft[i] = new TH1F(name,name,200,-0.1,1.9);
    sprintf(name,"thetamuAftRw[%d]",i);
    thetamuAftRw[i] = new TH1F(name,name,200,-0.1,1.9);

    sprintf(name,"q2qeB4[%d]",i);
    q2qeB4[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2qeB4Rw[%d]",i);
    q2qeB4Rw[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2qeAft[%d]",i);
    q2qeAft[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2qeAftRw[%d]",i);
    q2qeAftRw[i] = new TH1F(name,name,1200,0.0,40.0);

    sprintf(name,"q2B4[%d]",i);
    q2B4[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2B4Rw[%d]",i);
    q2B4Rw[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2Aft[%d]",i);
    q2Aft[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2AftRw[%d]",i);
    q2AftRw[i] = new TH1F(name,name,1200,0.0,40.0);

    sprintf(name,"xB4[%d]",i);
    xB4[i] = new TH1F(name,name,1200,0.0,30.0);
    sprintf(name,"xB4Rw[%d]",i);
    xB4Rw[i] = new TH1F(name,name,1200,0.0,30.0);
    sprintf(name,"xAft[%d]",i);
    xAft[i] = new TH1F(name,name,1200,0.0,30.0);
    sprintf(name,"xAftRw[%d]",i);
    xAftRw[i] = new TH1F(name,name,1200,0.0,30.0);

    sprintf(name,"yB4[%d]",i);
    yB4[i] = new TH1F(name,name,200,0.0,1.0);
    sprintf(name,"yB4Rw[%d]",i);
    yB4Rw[i] = new TH1F(name,name,200,0.0,1.0);
    sprintf(name,"yAft[%d]",i);
    yAft[i] = new TH1F(name,name,200,0.0,1.0);
    sprintf(name,"yAftRw[%d]",i);
    yAftRw[i] = new TH1F(name,name,200,0.0,1.0);

    sprintf(name,"w2B4[%d]",i);
    w2B4[i] = new TH1F(name,name,300,-5.0,10.0);
    sprintf(name,"w2B4Rw[%d]",i);
    w2B4Rw[i] = new TH1F(name,name,300,-5.0,10.0);
    sprintf(name,"w2Aft[%d]",i);
    w2Aft[i] = new TH1F(name,name,300,-5.0,10.0);
    sprintf(name,"w2AftRw[%d]",i);
    w2AftRw[i] = new TH1F(name,name,300,-5.0,10.0);

    // 2D (energies on x-axis where possible)

    sprintf(name,"thetamuVemuB4[%d]",i);
    thetamuVemuB4[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);
    sprintf(name,"thetamuVemuB4Rw[%d]",i);
    thetamuVemuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);
    sprintf(name,"thetamuVemuAft[%d]",i);
    thetamuVemuAft[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);
    sprintf(name,"thetamuVemuAftRw[%d]",i);
    thetamuVemuAftRw[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);

    sprintf(name,"thetamuVeshwB4[%d]",i);
    thetamuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);
    sprintf(name,"thetamuVeshwB4Rw[%d]",i);
    thetamuVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);
    sprintf(name,"thetamuVeshwAft[%d]",i);
    thetamuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);
    sprintf(name,"thetamuVeshwAftRw[%d]",i);
    thetamuVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);

    sprintf(name,"yVenuB4[%d]",i);
    yVenuB4[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    sprintf(name,"yVenuB4Rw[%d]",i);
    yVenuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    sprintf(name,"yVenuAft[%d]",i);
    yVenuAft[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    sprintf(name,"yVenuAftRw[%d]",i);
    yVenuAftRw[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    
    sprintf(name,"q2VenuB4[%d]",i);
    q2VenuB4[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);
    sprintf(name,"q2VenuB4Rw[%d]",i);
    q2VenuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);
    sprintf(name,"q2VenuAft[%d]",i);
    q2VenuAft[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);
    sprintf(name,"q2VenuAftRw[%d]",i);
    q2VenuAftRw[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);

    sprintf(name,"w2VenuB4[%d]",i);
    w2VenuB4[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);
    sprintf(name,"w2VenuB4Rw[%d]",i);
    w2VenuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);
    sprintf(name,"w2VenuAft[%d]",i);
    w2VenuAft[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);
    sprintf(name,"w2VenuAftRw[%d]",i);
    w2VenuAftRw[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);

    sprintf(name,"xVenuB4[%d]",i);
    xVenuB4[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);
    sprintf(name,"xVenuB4Rw[%d]",i);
    xVenuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);
    sprintf(name,"xVenuAft[%d]",i);
    xVenuAft[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);
    sprintf(name,"xVenuAftRw[%d]",i);
    xVenuAftRw[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);

    sprintf(name,"xVq2B4[%d]",i);
    xVq2B4[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    sprintf(name,"xVq2B4Rw[%d]",i);
    xVq2B4Rw[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    sprintf(name,"xVq2Aft[%d]",i);
    xVq2Aft[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    sprintf(name,"xVq2AftRw[%d]",i);
    xVq2AftRw[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    
    sprintf(name,"rmtotenuVenuB4[%d]",i);
    rmtotenuVenuB4[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotenuVenuB4Rw[%d]",i);
    rmtotenuVenuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotenuVenuAft[%d]",i);
    rmtotenuVenuAft[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotenuVenuAftRw[%d]",i);
    rmtotenuVenuAftRw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    
    sprintf(name,"rmtotemuVemuB4[%d]",i);
    rmtotemuVemuB4[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotemuVemuB4Rw[%d]",i);
    rmtotemuVemuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotemuVemuAft[%d]",i);
    rmtotemuVemuAft[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotemuVemuAftRw[%d]",i);
    rmtotemuVemuAftRw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);

    sprintf(name,"rmtoteshwVeshwB4[%d]",i);
    rmtoteshwVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,400,-10.0,10.0);
    sprintf(name,"rmtoteshwVeshwB4Rw[%d]",i);
    rmtoteshwVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,400,-10.0,10.0);
    sprintf(name,"rmtoteshwVeshwAft[%d]",i);
    rmtoteshwVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,400,-10.0,10.0);
    sprintf(name,"rmtoteshwVeshwAftRw[%d]",i);
    rmtoteshwVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,400,-10.0,10.0);
    
    sprintf(name,"rmtotthetamuVthetamuB4[%d]",i);
    rmtotthetamuVthetamuB4[i] = new TH2F(name,name,200,-0.1,1.9,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVthetamuB4Rw[%d]",i);
    rmtotthetamuVthetamuB4Rw[i] = new TH2F(name,name,200,-0.1,1.9,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVthetamuAft[%d]",i);
    rmtotthetamuVthetamuAft[i] = new TH2F(name,name,200,-0.1,1.9,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVthetamuAftRw[%d]",i);
    rmtotthetamuVthetamuAftRw[i] = new TH2F(name,name,200,-0.1,1.9,400,-10.0,10.0);
    
    sprintf(name,"rmtotthetamuVemuB4[%d]",i);
    rmtotthetamuVemuB4[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVemuB4Rw[%d]",i);
    rmtotthetamuVemuB4Rw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVemuAft[%d]",i);
    rmtotthetamuVemuAft[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
    sprintf(name,"rmtotthetamuVemuAftRw[%d]",i);
    rmtotthetamuVemuAftRw[i] = new TH2F(name,name,960,0.0,120.0,400,-10.0,10.0);
 
    sprintf(name,"enuVeshwB4[%d]",i);
    enuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuVeshwB4Rw[%d]",i);
    enuVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuVeshwAft[%d]",i);
    enuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuVeshwAftRw[%d]",i);
    enuVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"emuVeshwB4[%d]",i);
    emuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"emuVeshwB4Rw[%d]",i);
    emuVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"emuVeshwAft[%d]",i);
    emuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"emuVeshwAftRw[%d]",i);
    emuVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"enuqeVeshwB4[%d]",i);
    enuqeVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuqeVeshwB4Rw[%d]",i);
    enuqeVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuqeVeshwAft[%d]",i);
    enuqeVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuqeVeshwAftRw[%d]",i);
    enuqeVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"q2qeVeshwB4[%d]",i);
    q2qeVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2qeVeshwB4Rw[%d]",i);
    q2qeVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2qeVeshwAft[%d]",i);
    q2qeVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2qeVeshwAftRw[%d]",i);
    q2qeVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);

    sprintf(name,"q2VeshwB4[%d]",i);
    q2VeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2VeshwB4Rw[%d]",i);
    q2VeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2VeshwAft[%d]",i);
    q2VeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2VeshwAftRw[%d]",i);
    q2VeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);

    sprintf(name,"xVeshwB4[%d]",i);
    xVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);
    sprintf(name,"xVeshwB4Rw[%d]",i);
    xVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);
    sprintf(name,"xVeshwAft[%d]",i);
    xVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);
    sprintf(name,"xVeshwAftRw[%d]",i);
    xVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);

    sprintf(name,"yVeshwB4[%d]",i);
    yVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);
    sprintf(name,"yVeshwB4Rw[%d]",i);
    yVeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);
    sprintf(name,"yVeshwAft[%d]",i);
    yVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);
    sprintf(name,"yVeshwAftRw[%d]",i);
    yVeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);

    sprintf(name,"w2VeshwB4[%d]",i);
    w2VeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);
    sprintf(name,"w2VeshwB4Rw[%d]",i);
    w2VeshwB4Rw[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);
    sprintf(name,"w2VeshwAft[%d]",i);
    w2VeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);
    sprintf(name,"w2VeshwAftRw[%d]",i);
    w2VeshwAftRw[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);

  }
  
  ///////////////////////////////////////////////////

  dir->cd();

  Int_t nSnarl = 0;
  Bool_t fstSnarl = true;
  Int_t crntSnarl = 0;
  Int_t nRun = 0;
  Int_t crntRun = 0;

  for(Long64_t jentry=0;jentry<nentries;jentry++){
    Long64_t ientry = LoadTree(jentry);
    if(ientry<0) break;
    nb = fChain->GetEntry(jentry);
    nBytes += nb;

    // Run and snarl book keeping

    if(fstSnarl){
      crntSnarl = snarl;
      crntRun = run;
      nSnarl++;
      nRun++;
      fstSnarl = false;
    }
    else {
      if(crntRun!=run){
	crntSnarl = snarl;
	crntRun = run;
	nSnarl++;
	nRun++;
      }
      else if(crntSnarl!=snarl){
	crntSnarl = snarl;
	nSnarl++;
      }
    }

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

    // Make histogram files that correspond to specific processes

    if(processToKeep>=1001 && processToKeep<=1004){
      if(inu!=14) continue;
      if(cc_nc!=1) continue;
      if(process!=processToKeep) continue;
    }

    if(lookAtNC){
      if(inu!=14) continue;
      if(cc_nc!=0) continue;
    }

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

    // Good event cuts
      
    if(reco_enu>=120) continue;
    // Pitt ND fid vol etc.
    if(pass!=1) continue;
    // Track fit momentum error cut
    if(emu_meth==1 && fabs(trkeqp/trkqp)>0.3) continue;
    // Anti-neutrino cut
    if(trkqp>0) continue;
    // Rock muon cut
    if(reco_trk_vtxz<0.5) continue;
    // Runt event cut
    if((reco_enu - nvtp)/(reco_enu + nvtp)<0.95 && (nvtt*1e9)<50) continue;

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

    // Best approximation of 'standard' CC group good event cuts
      /*
    if(ntrack==0) continue; 
    if(reco_trk_vtxz<=1 || reco_trk_vtxz>=5) continue;
    if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))>=1) continue;
    if(pitt_trk_qual!=1) continue;
    if(reco_dircosneu<=0.6) continue;
    if(trkqp>0) continue;
      */

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

    // Do CC PID cut?
    if(selectCC){
      if(dave_cc_pid<-0.1) continue;
    }

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

    // Fiducial Sample Cuts
      
    // 'Anulus' 
    //if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))<0.4 || 
    //   (((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))>=100.0) continue;

    // Quadrants with beam spot 'chopped out'
    //if(reco_trk_vtxx>1.4885) continue; 
    //if(reco_trk_vtxy<=0.0) continue;  
    //if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))<0.04) continue; 

    // Z chunks
    if(reco_trk_vtxz<2.5 || reco_trk_vtxz>=3.5) continue;

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

    // Fill histograms before sample selection
    
    enuB4[0]->Fill(reco_enu);
    enuB4Rw[0]->Fill(reco_enu,neuweight*beamweight*hadweight);
    emuB4[0]->Fill(reco_emu);
    emuB4Rw[0]->Fill(reco_emu,neuweight*beamweight*hadweight);
    eshwB4[0]->Fill(reco_eshw);
    eshwB4Rw[0]->Fill(reco_eshw,neuweight*beamweight*hadweight);
    enuqeB4[0]->Fill(reco_qe_enu);
    enuqeB4Rw[0]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
    thetamuB4[0]->Fill(acos(reco_dircosneu));
    thetamuB4Rw[0]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
    q2qeB4[0]->Fill(reco_qe_q2);
    q2qeB4Rw[0]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
    q2B4[0]->Fill(reco_Q2);
    q2B4Rw[0]->Fill(reco_Q2,neuweight*beamweight*hadweight);
    xB4[0]->Fill(reco_x);
    xB4Rw[0]->Fill(reco_x,neuweight*beamweight*hadweight);
    yB4[0]->Fill(reco_y);
    yB4Rw[0]->Fill(reco_y,neuweight*beamweight*hadweight);
    w2B4[0]->Fill(reco_W2);
    w2B4Rw[0]->Fill(reco_W2,neuweight*beamweight*hadweight);

    thetamuVemuB4[0]->Fill(reco_emu,acos(reco_dircosneu));
    thetamuVemuB4Rw[0]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
    thetamuVeshwB4[0]->Fill(reco_eshw,acos(reco_dircosneu));
    thetamuVeshwB4Rw[0]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
    yVenuB4[0]->Fill(reco_enu,reco_y);
    yVenuB4Rw[0]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
    q2VenuB4[0]->Fill(reco_enu,reco_Q2);
    q2VenuB4Rw[0]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
    w2VenuB4[0]->Fill(reco_enu,reco_W2);
    w2VenuB4Rw[0]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
    xVenuB4[0]->Fill(reco_enu,reco_x);
    xVenuB4Rw[0]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
    xVq2B4[0]->Fill(log10(reco_x),log10(reco_Q2));
    xVq2B4Rw[0]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
    rmtotenuVenuB4[0]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
    rmtotenuVenuB4Rw[0]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
    rmtotemuVemuB4[0]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
    rmtotemuVemuB4Rw[0]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
    rmtoteshwVeshwB4[0]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
    rmtoteshwVeshwB4Rw[0]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
    rmtotthetamuVthetamuB4[0]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
    rmtotthetamuVthetamuB4Rw[0]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
    rmtotthetamuVemuB4[0]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
    rmtotthetamuVemuB4Rw[0]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
    
    trkVtxB4->Fill(reco_trk_vtxx,reco_trk_vtxy);

    enuVeshwB4[0]->Fill(reco_eshw,reco_enu);
    enuVeshwB4Rw[0]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
    emuVeshwB4[0]->Fill(reco_eshw,reco_emu);
    emuVeshwB4Rw[0]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
    enuqeVeshwB4[0]->Fill(reco_eshw,reco_qe_enu);
    enuqeVeshwB4Rw[0]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
    q2qeVeshwB4[0]->Fill(reco_eshw,reco_qe_q2);
    q2qeVeshwB4Rw[0]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
    q2VeshwB4[0]->Fill(reco_eshw,reco_Q2);
    q2VeshwB4Rw[0]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
    xVeshwB4[0]->Fill(reco_eshw,reco_x);
    xVeshwB4Rw[0]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
    yVeshwB4[0]->Fill(reco_eshw,reco_y);
    yVeshwB4Rw[0]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
    w2VeshwB4[0]->Fill(reco_eshw,reco_W2);
    w2VeshwB4Rw[0]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);

    if(reco_Q2>1 && reco_W2>4){
      enuB4[1]->Fill(reco_enu);
      enuB4Rw[1]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuB4[1]->Fill(reco_emu);
      emuB4Rw[1]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwB4[1]->Fill(reco_eshw);
      eshwB4Rw[1]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeB4[1]->Fill(reco_qe_enu);
      enuqeB4Rw[1]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuB4[1]->Fill(acos(reco_dircosneu));
      thetamuB4Rw[1]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeB4[1]->Fill(reco_qe_q2);
      q2qeB4Rw[1]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2B4[1]->Fill(reco_Q2);
      q2B4Rw[1]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xB4[1]->Fill(reco_x);
      xB4Rw[1]->Fill(reco_x,neuweight*beamweight*hadweight);
      yB4[1]->Fill(reco_y);
      yB4Rw[1]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2B4[1]->Fill(reco_W2);
      w2B4Rw[1]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuB4[1]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuB4Rw[1]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwB4[1]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwB4Rw[1]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuB4[1]->Fill(reco_enu,reco_y);
      yVenuB4Rw[1]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuB4[1]->Fill(reco_enu,reco_Q2);
      q2VenuB4Rw[1]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuB4[1]->Fill(reco_enu,reco_W2);
      w2VenuB4Rw[1]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuB4[1]->Fill(reco_enu,reco_x);
      xVenuB4Rw[1]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2B4[1]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2B4Rw[1]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuB4[1]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuB4Rw[1]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuB4[1]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuB4Rw[1]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwB4[1]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwB4Rw[1]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuB4[1]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuB4Rw[1]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuB4[1]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuB4Rw[1]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwB4[1]->Fill(reco_eshw,reco_enu);
      enuVeshwB4Rw[1]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwB4[1]->Fill(reco_eshw,reco_emu);
      emuVeshwB4Rw[1]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwB4[1]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwB4Rw[1]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwB4[1]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwB4Rw[1]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwB4[1]->Fill(reco_eshw,reco_Q2);
      q2VeshwB4Rw[1]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwB4[1]->Fill(reco_eshw,reco_x);
      xVeshwB4Rw[1]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwB4[1]->Fill(reco_eshw,reco_y);
      yVeshwB4Rw[1]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwB4[1]->Fill(reco_eshw,reco_W2);
      w2VeshwB4Rw[1]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_Q2<=1 && reco_W2>4){
      enuB4[2]->Fill(reco_enu);
      enuB4Rw[2]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuB4[2]->Fill(reco_emu);
      emuB4Rw[2]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwB4[2]->Fill(reco_eshw);
      eshwB4Rw[2]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeB4[2]->Fill(reco_qe_enu);
      enuqeB4Rw[2]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuB4[2]->Fill(acos(reco_dircosneu));
      thetamuB4Rw[2]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeB4[2]->Fill(reco_qe_q2);
      q2qeB4Rw[2]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2B4[2]->Fill(reco_Q2);
      q2B4Rw[2]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xB4[2]->Fill(reco_x);
      xB4Rw[2]->Fill(reco_x,neuweight*beamweight*hadweight);
      yB4[2]->Fill(reco_y);
      yB4Rw[2]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2B4[2]->Fill(reco_W2);
      w2B4Rw[2]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuB4[2]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuB4Rw[2]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwB4[2]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwB4Rw[2]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuB4[2]->Fill(reco_enu,reco_y);
      yVenuB4Rw[2]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuB4[2]->Fill(reco_enu,reco_Q2);
      q2VenuB4Rw[2]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuB4[2]->Fill(reco_enu,reco_W2);
      w2VenuB4Rw[2]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuB4[2]->Fill(reco_enu,reco_x);
      xVenuB4Rw[2]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2B4[2]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2B4Rw[2]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuB4[2]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuB4Rw[2]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuB4[2]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuB4Rw[2]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwB4[2]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwB4Rw[2]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuB4[2]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuB4Rw[2]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuB4[2]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuB4Rw[2]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwB4[2]->Fill(reco_eshw,reco_enu);
      enuVeshwB4Rw[2]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwB4[2]->Fill(reco_eshw,reco_emu);
      emuVeshwB4Rw[2]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwB4[2]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwB4Rw[2]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwB4[2]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwB4Rw[2]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwB4[2]->Fill(reco_eshw,reco_Q2);
      q2VeshwB4Rw[2]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwB4[2]->Fill(reco_eshw,reco_x);
      xVeshwB4Rw[2]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwB4[2]->Fill(reco_eshw,reco_y);
      yVeshwB4Rw[2]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwB4[2]->Fill(reco_eshw,reco_W2);
      w2VeshwB4Rw[2]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_W2<1.69){
      enuB4[3]->Fill(reco_enu);
      enuB4Rw[3]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuB4[3]->Fill(reco_emu);
      emuB4Rw[3]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwB4[3]->Fill(reco_eshw);
      eshwB4Rw[3]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeB4[3]->Fill(reco_qe_enu);
      enuqeB4Rw[3]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuB4[3]->Fill(acos(reco_dircosneu));
      thetamuB4Rw[3]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeB4[3]->Fill(reco_qe_q2);
      q2qeB4Rw[3]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2B4[3]->Fill(reco_Q2);
      q2B4Rw[3]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xB4[3]->Fill(reco_x);
      xB4Rw[3]->Fill(reco_x,neuweight*beamweight*hadweight);
      yB4[3]->Fill(reco_y);
      yB4Rw[3]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2B4[3]->Fill(reco_W2);
      w2B4Rw[3]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuB4[3]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuB4Rw[3]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwB4[3]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwB4Rw[3]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuB4[3]->Fill(reco_enu,reco_y);
      yVenuB4Rw[3]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuB4[3]->Fill(reco_enu,reco_Q2);
      q2VenuB4Rw[3]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuB4[3]->Fill(reco_enu,reco_W2);
      w2VenuB4Rw[3]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuB4[3]->Fill(reco_enu,reco_x);
      xVenuB4Rw[3]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2B4[3]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2B4Rw[3]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuB4[3]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuB4Rw[3]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuB4[3]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuB4Rw[3]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwB4[3]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwB4Rw[3]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuB4[3]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuB4Rw[3]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuB4[3]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuB4Rw[3]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwB4[3]->Fill(reco_eshw,reco_enu);
      enuVeshwB4Rw[3]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwB4[3]->Fill(reco_eshw,reco_emu);
      emuVeshwB4Rw[3]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwB4[3]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwB4Rw[3]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwB4[3]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwB4Rw[3]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwB4[3]->Fill(reco_eshw,reco_Q2);
      q2VeshwB4Rw[3]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwB4[3]->Fill(reco_eshw,reco_x);
      xVeshwB4Rw[3]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwB4[3]->Fill(reco_eshw,reco_y);
      yVeshwB4Rw[3]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwB4[3]->Fill(reco_eshw,reco_W2);
      w2VeshwB4Rw[3]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_W2>=1.69 && reco_W2<=4){
      enuB4[4]->Fill(reco_enu);
      enuB4Rw[4]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuB4[4]->Fill(reco_emu);
      emuB4Rw[4]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwB4[4]->Fill(reco_eshw);
      eshwB4Rw[4]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeB4[4]->Fill(reco_qe_enu);
      enuqeB4Rw[4]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuB4[4]->Fill(acos(reco_dircosneu));
      thetamuB4Rw[4]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeB4[4]->Fill(reco_qe_q2);
      q2qeB4Rw[4]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2B4[4]->Fill(reco_Q2);
      q2B4Rw[4]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xB4[4]->Fill(reco_x);
      xB4Rw[4]->Fill(reco_x,neuweight*beamweight*hadweight);
      yB4[4]->Fill(reco_y);
      yB4Rw[4]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2B4[4]->Fill(reco_W2);
      w2B4Rw[4]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuB4[4]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuB4Rw[4]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwB4[4]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwB4Rw[4]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuB4[4]->Fill(reco_enu,reco_y);
      yVenuB4Rw[4]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuB4[4]->Fill(reco_enu,reco_Q2);
      q2VenuB4Rw[4]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuB4[4]->Fill(reco_enu,reco_W2);
      w2VenuB4Rw[4]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuB4[4]->Fill(reco_enu,reco_x);
      xVenuB4Rw[4]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2B4[4]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2B4Rw[4]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuB4[4]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuB4Rw[4]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuB4[4]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuB4Rw[4]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwB4[4]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwB4Rw[4]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuB4[4]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuB4Rw[4]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuB4[4]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuB4Rw[4]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwB4[4]->Fill(reco_eshw,reco_enu);
      enuVeshwB4Rw[4]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwB4[4]->Fill(reco_eshw,reco_emu);
      emuVeshwB4Rw[4]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwB4[4]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwB4Rw[4]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwB4[4]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwB4Rw[4]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwB4[4]->Fill(reco_eshw,reco_Q2);
      q2VeshwB4Rw[4]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwB4[4]->Fill(reco_eshw,reco_x);
      xVeshwB4Rw[4]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwB4[4]->Fill(reco_eshw,reco_y);
      yVeshwB4Rw[4]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwB4[4]->Fill(reco_eshw,reco_W2);
      w2VeshwB4Rw[4]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }

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

    // Select sample

    if(reco_eshw>0.0) continue;

    //if(med_qe_pass!=1) continue;

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

    // Fill histograms after sample selection

    enuAft[0]->Fill(reco_enu);
    enuAftRw[0]->Fill(reco_enu,neuweight*beamweight*hadweight);
    emuAft[0]->Fill(reco_emu);
    emuAftRw[0]->Fill(reco_emu,neuweight*beamweight*hadweight);
    eshwAft[0]->Fill(reco_eshw);
    eshwAftRw[0]->Fill(reco_eshw,neuweight*beamweight*hadweight);
    enuqeAft[0]->Fill(reco_qe_enu);
    enuqeAftRw[0]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
    thetamuAft[0]->Fill(acos(reco_dircosneu));
    thetamuAftRw[0]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
    q2qeAft[0]->Fill(reco_qe_q2);
    q2qeAftRw[0]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
    q2Aft[0]->Fill(reco_Q2);
    q2AftRw[0]->Fill(reco_Q2,neuweight*beamweight*hadweight);
    xAft[0]->Fill(reco_x);
    xAftRw[0]->Fill(reco_x,neuweight*beamweight*hadweight);
    yAft[0]->Fill(reco_y);
    yAftRw[0]->Fill(reco_y,neuweight*beamweight*hadweight);
    w2Aft[0]->Fill(reco_W2);
    w2AftRw[0]->Fill(reco_W2,neuweight*beamweight*hadweight);

    thetamuVemuAft[0]->Fill(reco_emu,acos(reco_dircosneu));
    thetamuVemuAftRw[0]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
    thetamuVeshwAft[0]->Fill(reco_eshw,acos(reco_dircosneu));
    thetamuVeshwAftRw[0]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
    yVenuAft[0]->Fill(reco_enu,reco_y);
    yVenuAftRw[0]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
    q2VenuAft[0]->Fill(reco_enu,reco_Q2);
    q2VenuAftRw[0]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
    w2VenuAft[0]->Fill(reco_enu,reco_W2);
    w2VenuAftRw[0]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
    xVenuAft[0]->Fill(reco_enu,reco_x);
    xVenuAftRw[0]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
    xVq2Aft[0]->Fill(log10(reco_x),log10(reco_Q2));
    xVq2AftRw[0]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
    rmtotenuVenuAft[0]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
    rmtotenuVenuAftRw[0]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
    rmtotemuVemuAft[0]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
    rmtotemuVemuAftRw[0]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
    rmtoteshwVeshwAft[0]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
    rmtoteshwVeshwAftRw[0]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
    rmtotthetamuVthetamuAft[0]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
    rmtotthetamuVthetamuAftRw[0]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
    rmtotthetamuVemuAft[0]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
    rmtotthetamuVemuAftRw[0]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
    
    trkVtxAft->Fill(reco_trk_vtxx,reco_trk_vtxy);

    enuVeshwAft[0]->Fill(reco_eshw,reco_enu);
    enuVeshwAftRw[0]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
    emuVeshwAft[0]->Fill(reco_eshw,reco_emu);
    emuVeshwAftRw[0]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
    enuqeVeshwAft[0]->Fill(reco_eshw,reco_qe_enu);
    enuqeVeshwAftRw[0]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
    q2qeVeshwAft[0]->Fill(reco_eshw,reco_qe_q2);
    q2qeVeshwAftRw[0]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
    q2VeshwAft[0]->Fill(reco_eshw,reco_Q2);
    q2VeshwAftRw[0]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
    xVeshwAft[0]->Fill(reco_eshw,reco_x);
    xVeshwAftRw[0]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
    yVeshwAft[0]->Fill(reco_eshw,reco_y);
    yVeshwAftRw[0]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
    w2VeshwAft[0]->Fill(reco_eshw,reco_W2);
    w2VeshwAftRw[0]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);

    if(reco_Q2>1 && reco_W2>4){
      enuAft[1]->Fill(reco_enu);
      enuAftRw[1]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuAft[1]->Fill(reco_emu);
      emuAftRw[1]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwAft[1]->Fill(reco_eshw);
      eshwAftRw[1]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeAft[1]->Fill(reco_qe_enu);
      enuqeAftRw[1]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuAft[1]->Fill(acos(reco_dircosneu));
      thetamuAftRw[1]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeAft[1]->Fill(reco_qe_q2);
      q2qeAftRw[1]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2Aft[1]->Fill(reco_Q2);
      q2AftRw[1]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xAft[1]->Fill(reco_x);
      xAftRw[1]->Fill(reco_x,neuweight*beamweight*hadweight);
      yAft[1]->Fill(reco_y);
      yAftRw[1]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2Aft[1]->Fill(reco_W2);
      w2AftRw[1]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuAft[1]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuAftRw[1]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwAft[1]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwAftRw[1]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuAft[1]->Fill(reco_enu,reco_y);
      yVenuAftRw[1]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuAft[1]->Fill(reco_enu,reco_Q2);
      q2VenuAftRw[1]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuAft[1]->Fill(reco_enu,reco_W2);
      w2VenuAftRw[1]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuAft[1]->Fill(reco_enu,reco_x);
      xVenuAftRw[1]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2Aft[1]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2AftRw[1]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuAft[1]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuAftRw[1]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuAft[1]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuAftRw[1]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwAft[1]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwAftRw[1]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuAft[1]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuAftRw[1]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuAft[1]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuAftRw[1]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwAft[1]->Fill(reco_eshw,reco_enu);
      enuVeshwAftRw[1]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwAft[1]->Fill(reco_eshw,reco_emu);
      emuVeshwAftRw[1]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwAft[1]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwAftRw[1]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwAft[1]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwAftRw[1]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwAft[1]->Fill(reco_eshw,reco_Q2);
      q2VeshwAftRw[1]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwAft[1]->Fill(reco_eshw,reco_x);
      xVeshwAftRw[1]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwAft[1]->Fill(reco_eshw,reco_y);
      yVeshwAftRw[1]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwAft[1]->Fill(reco_eshw,reco_W2);
      w2VeshwAftRw[1]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_Q2<=1 && reco_W2>4){
      enuAft[2]->Fill(reco_enu);
      enuAftRw[2]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuAft[2]->Fill(reco_emu);
      emuAftRw[2]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwAft[2]->Fill(reco_eshw);
      eshwAftRw[2]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeAft[2]->Fill(reco_qe_enu);
      enuqeAftRw[2]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuAft[2]->Fill(acos(reco_dircosneu));
      thetamuAftRw[2]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeAft[2]->Fill(reco_qe_q2);
      q2qeAftRw[2]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2Aft[2]->Fill(reco_Q2);
      q2AftRw[2]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xAft[2]->Fill(reco_x);
      xAftRw[2]->Fill(reco_x,neuweight*beamweight*hadweight);
      yAft[2]->Fill(reco_y);
      yAftRw[2]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2Aft[2]->Fill(reco_W2);
      w2AftRw[2]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuAft[2]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuAftRw[2]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwAft[2]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwAftRw[2]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuAft[2]->Fill(reco_enu,reco_y);
      yVenuAftRw[2]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuAft[2]->Fill(reco_enu,reco_Q2);
      q2VenuAftRw[2]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuAft[2]->Fill(reco_enu,reco_W2);
      w2VenuAftRw[2]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuAft[2]->Fill(reco_enu,reco_x);
      xVenuAftRw[2]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2Aft[2]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2AftRw[2]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuAft[2]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuAftRw[2]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuAft[2]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuAftRw[2]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwAft[2]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwAftRw[2]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuAft[2]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuAftRw[2]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuAft[2]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuAftRw[2]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwAft[2]->Fill(reco_eshw,reco_enu);
      enuVeshwAftRw[2]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwAft[2]->Fill(reco_eshw,reco_emu);
      emuVeshwAftRw[2]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwAft[2]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwAftRw[2]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwAft[2]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwAftRw[2]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwAft[2]->Fill(reco_eshw,reco_Q2);
      q2VeshwAftRw[2]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwAft[2]->Fill(reco_eshw,reco_x);
      xVeshwAftRw[2]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwAft[2]->Fill(reco_eshw,reco_y);
      yVeshwAftRw[2]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwAft[2]->Fill(reco_eshw,reco_W2);
      w2VeshwAftRw[2]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_W2<1.69){
      enuAft[3]->Fill(reco_enu);
      enuAftRw[3]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuAft[3]->Fill(reco_emu);
      emuAftRw[3]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwAft[3]->Fill(reco_eshw);
      eshwAftRw[3]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeAft[3]->Fill(reco_qe_enu);
      enuqeAftRw[3]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuAft[3]->Fill(acos(reco_dircosneu));
      thetamuAftRw[3]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeAft[3]->Fill(reco_qe_q2);
      q2qeAftRw[3]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2Aft[3]->Fill(reco_Q2);
      q2AftRw[3]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xAft[3]->Fill(reco_x);
      xAftRw[3]->Fill(reco_x,neuweight*beamweight*hadweight);
      yAft[3]->Fill(reco_y);
      yAftRw[3]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2Aft[3]->Fill(reco_W2);
      w2AftRw[3]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuAft[3]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuAftRw[3]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwAft[3]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwAftRw[3]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuAft[3]->Fill(reco_enu,reco_y);
      yVenuAftRw[3]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuAft[3]->Fill(reco_enu,reco_Q2);
      q2VenuAftRw[3]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuAft[3]->Fill(reco_enu,reco_W2);
      w2VenuAftRw[3]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuAft[3]->Fill(reco_enu,reco_x);
      xVenuAftRw[3]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2Aft[3]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2AftRw[3]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuAft[3]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuAftRw[3]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuAft[3]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuAftRw[3]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwAft[3]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwAftRw[3]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuAft[3]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuAftRw[3]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuAft[3]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuAftRw[3]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwAft[3]->Fill(reco_eshw,reco_enu);
      enuVeshwAftRw[3]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwAft[3]->Fill(reco_eshw,reco_emu);
      emuVeshwAftRw[3]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwAft[3]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwAftRw[3]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwAft[3]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwAftRw[3]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwAft[3]->Fill(reco_eshw,reco_Q2);
      q2VeshwAftRw[3]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwAft[3]->Fill(reco_eshw,reco_x);
      xVeshwAftRw[3]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwAft[3]->Fill(reco_eshw,reco_y);
      yVeshwAftRw[3]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwAft[3]->Fill(reco_eshw,reco_W2);
      w2VeshwAftRw[3]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }
    if(reco_W2>=1.69 && reco_W2<=4){
      enuAft[4]->Fill(reco_enu);
      enuAftRw[4]->Fill(reco_enu,neuweight*beamweight*hadweight);
      emuAft[4]->Fill(reco_emu);
      emuAftRw[4]->Fill(reco_emu,neuweight*beamweight*hadweight);
      eshwAft[4]->Fill(reco_eshw);
      eshwAftRw[4]->Fill(reco_eshw,neuweight*beamweight*hadweight);
      enuqeAft[4]->Fill(reco_qe_enu);
      enuqeAftRw[4]->Fill(reco_qe_enu,neuweight*beamweight*hadweight);
      thetamuAft[4]->Fill(acos(reco_dircosneu));
      thetamuAftRw[4]->Fill(acos(reco_dircosneu),neuweight*beamweight*hadweight);
      q2qeAft[4]->Fill(reco_qe_q2);
      q2qeAftRw[4]->Fill(reco_qe_q2,neuweight*beamweight*hadweight);
      q2Aft[4]->Fill(reco_Q2);
      q2AftRw[4]->Fill(reco_Q2,neuweight*beamweight*hadweight);
      xAft[4]->Fill(reco_x);
      xAftRw[4]->Fill(reco_x,neuweight*beamweight*hadweight);
      yAft[4]->Fill(reco_y);
      yAftRw[4]->Fill(reco_y,neuweight*beamweight*hadweight);
      w2Aft[4]->Fill(reco_W2);
      w2AftRw[4]->Fill(reco_W2,neuweight*beamweight*hadweight);
      
      thetamuVemuAft[4]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVemuAftRw[4]->Fill(reco_emu,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      thetamuVeshwAft[4]->Fill(reco_eshw,acos(reco_dircosneu));
      thetamuVeshwAftRw[4]->Fill(reco_eshw,acos(reco_dircosneu),neuweight*beamweight*hadweight);
      yVenuAft[4]->Fill(reco_enu,reco_y);
      yVenuAftRw[4]->Fill(reco_enu,reco_y,neuweight*beamweight*hadweight);
      q2VenuAft[4]->Fill(reco_enu,reco_Q2);
      q2VenuAftRw[4]->Fill(reco_enu,reco_Q2,neuweight*beamweight*hadweight);
      w2VenuAft[4]->Fill(reco_enu,reco_W2);
      w2VenuAftRw[4]->Fill(reco_enu,reco_W2,neuweight*beamweight*hadweight);
      xVenuAft[4]->Fill(reco_enu,reco_x);
      xVenuAftRw[4]->Fill(reco_enu,reco_x,neuweight*beamweight*hadweight);
      xVq2Aft[4]->Fill(log10(reco_x),log10(reco_Q2));
      xVq2AftRw[4]->Fill(log10(reco_x),log10(reco_Q2),neuweight*beamweight*hadweight);
      rmtotenuVenuAft[4]->Fill(true_enu,((reco_enu-true_enu)/true_enu));
      rmtotenuVenuAftRw[4]->Fill(true_enu,((reco_enu-true_enu)/true_enu),neuweight*beamweight*hadweight);
      rmtotemuVemuAft[4]->Fill(true_emu,((reco_emu-true_emu)/true_emu));
      rmtotemuVemuAftRw[4]->Fill(true_emu,((reco_emu-true_emu)/true_emu),neuweight*beamweight*hadweight);
      rmtoteshwVeshwAft[4]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw));
      rmtoteshwVeshwAftRw[4]->Fill(true_eshw,((reco_eshw-true_eshw)/true_eshw),neuweight*beamweight*hadweight);
      rmtotthetamuVthetamuAft[4]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVthetamuAftRw[4]->Fill(acos(true_dircosneu),((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);
      rmtotthetamuVemuAft[4]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)));
      rmtotthetamuVemuAftRw[4]->Fill(true_emu,((acos(reco_dircosneu)-acos(true_dircosneu))/acos(true_dircosneu)),neuweight*beamweight*hadweight);

      enuVeshwAft[4]->Fill(reco_eshw,reco_enu);
      enuVeshwAftRw[4]->Fill(reco_eshw,reco_enu,neuweight*beamweight*hadweight);
      emuVeshwAft[4]->Fill(reco_eshw,reco_emu);
      emuVeshwAftRw[4]->Fill(reco_eshw,reco_emu,neuweight*beamweight*hadweight);
      enuqeVeshwAft[4]->Fill(reco_eshw,reco_qe_enu);
      enuqeVeshwAftRw[4]->Fill(reco_eshw,reco_qe_enu,neuweight*beamweight*hadweight);
      q2qeVeshwAft[4]->Fill(reco_eshw,reco_qe_q2);
      q2qeVeshwAftRw[4]->Fill(reco_eshw,reco_qe_q2,neuweight*beamweight*hadweight);
      q2VeshwAft[4]->Fill(reco_eshw,reco_Q2);
      q2VeshwAftRw[4]->Fill(reco_eshw,reco_Q2,neuweight*beamweight*hadweight);
      xVeshwAft[4]->Fill(reco_eshw,reco_x);
      xVeshwAftRw[4]->Fill(reco_eshw,reco_x,neuweight*beamweight*hadweight);
      yVeshwAft[4]->Fill(reco_eshw,reco_y);
      yVeshwAftRw[4]->Fill(reco_eshw,reco_y,neuweight*beamweight*hadweight);
      w2VeshwAft[4]->Fill(reco_eshw,reco_W2);
      w2VeshwAftRw[4]->Fill(reco_eshw,reco_W2,neuweight*beamweight*hadweight);
    }

  }

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

  // Output histograms to file

  for(Int_t i=0;i<5;i++){
      
    enuB4[i]->SetDirectory(fout);
    enuB4Rw[i]->SetDirectory(fout);
    enuAft[i]->SetDirectory(fout);
    enuAftRw[i]->SetDirectory(fout);
    
    emuB4[i]->SetDirectory(fout);
    emuB4Rw[i]->SetDirectory(fout);
    emuAft[i]->SetDirectory(fout);
    emuAftRw[i]->SetDirectory(fout);
      
    eshwB4[i]->SetDirectory(fout);
    eshwB4Rw[i]->SetDirectory(fout);
    eshwAft[i]->SetDirectory(fout);
    eshwAftRw[i]->SetDirectory(fout);
      
    enuqeB4[i]->SetDirectory(fout);
    enuqeB4Rw[i]->SetDirectory(fout);
    enuqeAft[i]->SetDirectory(fout);
    enuqeAftRw[i]->SetDirectory(fout);
      
    thetamuB4[i]->SetDirectory(fout);
    thetamuB4Rw[i]->SetDirectory(fout);
    thetamuAft[i]->SetDirectory(fout);
    thetamuAftRw[i]->SetDirectory(fout);
      
    q2qeB4[i]->SetDirectory(fout);
    q2qeB4Rw[i]->SetDirectory(fout);
    q2qeAft[i]->SetDirectory(fout);
    q2qeAftRw[i]->SetDirectory(fout);
      
    q2B4[i]->SetDirectory(fout);
    q2B4Rw[i]->SetDirectory(fout);
    q2Aft[i]->SetDirectory(fout);
    q2AftRw[i]->SetDirectory(fout);
      
    xB4[i]->SetDirectory(fout);
    xB4Rw[i]->SetDirectory(fout);
    xAft[i]->SetDirectory(fout);
    xAftRw[i]->SetDirectory(fout);
      
    yB4[i]->SetDirectory(fout);
    yB4Rw[i]->SetDirectory(fout);
    yAft[i]->SetDirectory(fout);
    yAftRw[i]->SetDirectory(fout);
      
    w2B4[i]->SetDirectory(fout);
    w2B4Rw[i]->SetDirectory(fout);
    w2Aft[i]->SetDirectory(fout);
    w2AftRw[i]->SetDirectory(fout);

    thetamuVemuB4[i]->SetDirectory(fout);
    thetamuVemuB4Rw[i]->SetDirectory(fout);
    thetamuVemuAft[i]->SetDirectory(fout);
    thetamuVemuAftRw[i]->SetDirectory(fout);
      
    thetamuVeshwB4[i]->SetDirectory(fout);
    thetamuVeshwB4Rw[i]->SetDirectory(fout);
    thetamuVeshwAft[i]->SetDirectory(fout);
    thetamuVeshwAftRw[i]->SetDirectory(fout);
      
    yVenuB4[i]->SetDirectory(fout);
    yVenuB4Rw[i]->SetDirectory(fout);
    yVenuAft[i]->SetDirectory(fout);
    yVenuAftRw[i]->SetDirectory(fout);
      
    q2VenuB4[i]->SetDirectory(fout);
    q2VenuB4Rw[i]->SetDirectory(fout);
    q2VenuAft[i]->SetDirectory(fout);
    q2VenuAftRw[i]->SetDirectory(fout);

    w2VenuB4[i]->SetDirectory(fout);
    w2VenuB4Rw[i]->SetDirectory(fout);
    w2VenuAft[i]->SetDirectory(fout);
    w2VenuAftRw[i]->SetDirectory(fout);

    xVenuB4[i]->SetDirectory(fout);
    xVenuB4Rw[i]->SetDirectory(fout);
    xVenuAft[i]->SetDirectory(fout);
    xVenuAftRw[i]->SetDirectory(fout);
      
    xVq2B4[i]->SetDirectory(fout);
    xVq2B4Rw[i]->SetDirectory(fout);
    xVq2Aft[i]->SetDirectory(fout);
    xVq2AftRw[i]->SetDirectory(fout);
      
    rmtotenuVenuB4[i]->SetDirectory(fout);
    rmtotenuVenuB4Rw[i]->SetDirectory(fout);
    rmtotenuVenuAft[i]->SetDirectory(fout);
    rmtotenuVenuAftRw[i]->SetDirectory(fout);
      
    rmtoteshwVeshwB4[i]->SetDirectory(fout);
    rmtoteshwVeshwB4Rw[i]->SetDirectory(fout);
    rmtoteshwVeshwAft[i]->SetDirectory(fout);
    rmtoteshwVeshwAftRw[i]->SetDirectory(fout);
      
    rmtotemuVemuB4[i]->SetDirectory(fout);
    rmtotemuVemuB4Rw[i]->SetDirectory(fout);
    rmtotemuVemuAft[i]->SetDirectory(fout);
    rmtotemuVemuAftRw[i]->SetDirectory(fout);
      
    rmtotthetamuVthetamuB4[i]->SetDirectory(fout);
    rmtotthetamuVthetamuB4Rw[i]->SetDirectory(fout);
    rmtotthetamuVthetamuAft[i]->SetDirectory(fout);
    rmtotthetamuVthetamuAftRw[i]->SetDirectory(fout);
      
    rmtotthetamuVemuB4[i]->SetDirectory(fout);
    rmtotthetamuVemuB4Rw[i]->SetDirectory(fout);
    rmtotthetamuVemuAft[i]->SetDirectory(fout);
    rmtotthetamuVemuAftRw[i]->SetDirectory(fout);

    enuVeshwB4[i]->SetDirectory(fout);
    enuVeshwB4Rw[i]->SetDirectory(fout);
    enuVeshwAft[i]->SetDirectory(fout);
    enuVeshwAftRw[i]->SetDirectory(fout);
    
    emuVeshwB4[i]->SetDirectory(fout);
    emuVeshwB4Rw[i]->SetDirectory(fout);
    emuVeshwAft[i]->SetDirectory(fout);
    emuVeshwAftRw[i]->SetDirectory(fout);
      
    enuqeVeshwB4[i]->SetDirectory(fout);
    enuqeVeshwB4Rw[i]->SetDirectory(fout);
    enuqeVeshwAft[i]->SetDirectory(fout);
    enuqeVeshwAftRw[i]->SetDirectory(fout);
      
    q2qeVeshwB4[i]->SetDirectory(fout);
    q2qeVeshwB4Rw[i]->SetDirectory(fout);
    q2qeVeshwAft[i]->SetDirectory(fout);
    q2qeVeshwAftRw[i]->SetDirectory(fout);
      
    q2VeshwB4[i]->SetDirectory(fout);
    q2VeshwB4Rw[i]->SetDirectory(fout);
    q2VeshwAft[i]->SetDirectory(fout);
    q2VeshwAftRw[i]->SetDirectory(fout);
      
    xVeshwB4[i]->SetDirectory(fout);
    xVeshwB4Rw[i]->SetDirectory(fout);
    xVeshwAft[i]->SetDirectory(fout);
    xVeshwAftRw[i]->SetDirectory(fout);
      
    yVeshwB4[i]->SetDirectory(fout);
    yVeshwB4Rw[i]->SetDirectory(fout);
    yVeshwAft[i]->SetDirectory(fout);
    yVeshwAftRw[i]->SetDirectory(fout);
      
    w2VeshwB4[i]->SetDirectory(fout);
    w2VeshwB4Rw[i]->SetDirectory(fout);
    w2VeshwAft[i]->SetDirectory(fout);
    w2VeshwAftRw[i]->SetDirectory(fout);

  }
  
  trkVtxB4->SetDirectory(fout);
  trkVtxAft->SetDirectory(fout);

  fout->Write();
  
  cout << "Number of snarls:   " << nSnarl << endl;
  cout << "Number of runs:     " << nRun << endl;

}

//***************************************************
void MedPanAna::AnalyseData(Int_t beamConfig, Bool_t selectCC, const char* tag)
{

  if(fChain==0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  Int_t nBytes = 0, nb = 0;

  TFile* fout = new TFile(tag,"recreate");
  TDirectory* dir = gDirectory;

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

  // Declare histograms

  // Histogram array convention:
  //   [0] - all events
  //   [1] - 'safe' Dis    Q2>1,  W2>4
  //   [2] - 'low q2' Dis  Q2<=1, W2>4
  //   [3] - QE/Delta      W2<1.69
  //   [4] - 'transition'  W2>=1.69 && W2<=4

  // 1D

  TH1F* enuB4[5];
  TH1F* enuAft[5];

  TH1F* emuB4[5];
  TH1F* emuAft[5];

  TH1F* eshwB4[5];
  TH1F* eshwAft[5];

  TH1F* enuqeB4[5];
  TH1F* enuqeAft[5];

  TH1F* thetamuB4[5];
  TH1F* thetamuAft[5];

  TH1F* q2qeB4[5];
  TH1F* q2qeAft[5];

  TH1F* q2B4[5];
  TH1F* q2Aft[5];

  TH1F* xB4[5];
  TH1F* xAft[5];

  TH1F* yB4[5];
  TH1F* yAft[5];

  TH1F* w2B4[5];
  TH1F* w2Aft[5];

  // 2D

  TH2F* thetamuVemuB4[5];
  TH2F* thetamuVemuAft[5];

  TH2F* thetamuVeshwB4[5];
  TH2F* thetamuVeshwAft[5];

  TH2F* yVenuB4[5];
  TH2F* yVenuAft[5];

  TH2F* q2VenuB4[5];
  TH2F* q2VenuAft[5];

  TH2F* w2VenuB4[5];
  TH2F* w2VenuAft[5];

  TH2F* xVenuB4[5];
  TH2F* xVenuAft[5];

  TH2F* xVq2B4[5];
  TH2F* xVq2Aft[5];

  TH2F* enuVeshwB4[5];
  TH2F* enuVeshwAft[5];

  TH2F* emuVeshwB4[5];
  TH2F* emuVeshwAft[5];

  TH2F* enuqeVeshwB4[5];
  TH2F* enuqeVeshwAft[5];

  TH2F* q2qeVeshwB4[5];
  TH2F* q2qeVeshwAft[5];

  TH2F* q2VeshwB4[5];
  TH2F* q2VeshwAft[5];

  TH2F* xVeshwB4[5];
  TH2F* xVeshwAft[5];

  TH2F* yVeshwB4[5];
  TH2F* yVeshwAft[5];

  TH2F* w2VeshwB4[5];
  TH2F* w2VeshwAft[5];

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

  TH2F* trkVtxB4 = new TH2F("trkVtxB4","trkVtxB4",800,0.5,2.5,600,-1.5,1.5);
  TH2F* trkVtxAft = new TH2F("trkVtxAft","trkVtxAft",800,0.5,2.5,600,-1.5,1.5);

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

  // Initialise histograms

  char name[50];
  for(Int_t i=0;i<5;i++){
    
    // 1D
    
    sprintf(name,"enuB4[%d]",i);
    enuB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuAft[%d]",i);
    enuAft[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"emuB4[%d]",i);
    emuB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"emuAft[%d]",i);
    emuAft[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"eshwB4[%d]",i);
    eshwB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"eshwAft[%d]",i);
    eshwAft[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"enuqeB4[%d]",i);
    enuqeB4[i] = new TH1F(name,name,960,0,120);
    sprintf(name,"enuqeAft[%d]",i);
    enuqeAft[i] = new TH1F(name,name,960,0,120);

    sprintf(name,"thetamuB4[%d]",i);
    thetamuB4[i] = new TH1F(name,name,200,-0.1,1.9);
    sprintf(name,"thetamuAft[%d]",i);
    thetamuAft[i] = new TH1F(name,name,200,-0.1,1.9);

    sprintf(name,"q2qeB4[%d]",i);
    q2qeB4[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2qeAft[%d]",i);
    q2qeAft[i] = new TH1F(name,name,1200,0.0,40.0);

    sprintf(name,"q2B4[%d]",i);
    q2B4[i] = new TH1F(name,name,1200,0.0,40.0);
    sprintf(name,"q2Aft[%d]",i);
    q2Aft[i] = new TH1F(name,name,1200,0.0,40.0);

    sprintf(name,"xB4[%d]",i);
    xB4[i] = new TH1F(name,name,1200,0.0,30.0);
    sprintf(name,"xAft[%d]",i);
    xAft[i] = new TH1F(name,name,1200,0.0,30.0);

    sprintf(name,"yB4[%d]",i);
    yB4[i] = new TH1F(name,name,200,0.0,1.0);
    sprintf(name,"yAft[%d]",i);
    yAft[i] = new TH1F(name,name,200,0.0,1.0);

    sprintf(name,"w2B4[%d]",i);
    w2B4[i] = new TH1F(name,name,300,-5.0,10.0);
    sprintf(name,"w2Aft[%d]",i);
    w2Aft[i] = new TH1F(name,name,300,-5.0,10.0);

    // 2D (energies on x-axis where possible)

    sprintf(name,"thetamuVemuB4[%d]",i);
    thetamuVemuB4[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);
    sprintf(name,"thetamuVemuAft[%d]",i);
    thetamuVemuAft[i] = new TH2F(name,name,960,0.0,120.0,200,-0.1,1.9);

    sprintf(name,"thetamuVeshwB4[%d]",i);
    thetamuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);
    sprintf(name,"thetamuVeshwAft[%d]",i);
    thetamuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,200,-0.1,1.9);

    sprintf(name,"yVenuB4[%d]",i);
    yVenuB4[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    sprintf(name,"yVenuAft[%d]",i);
    yVenuAft[i] = new TH2F(name,name,960,0.0,120.0,200,0.0,1.0);
    
    sprintf(name,"q2VenuB4[%d]",i);
    q2VenuB4[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);
    sprintf(name,"q2VenuAft[%d]",i);
    q2VenuAft[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,40.0);

    sprintf(name,"w2VenuB4[%d]",i);
    w2VenuB4[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);
    sprintf(name,"w2VenuAft[%d]",i);
    w2VenuAft[i] = new TH2F(name,name,960,0.0,120.0,300,-5.0,10.0);

    sprintf(name,"xVenuB4[%d]",i);
    xVenuB4[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);
    sprintf(name,"xVenuAft[%d]",i);
    xVenuAft[i] = new TH2F(name,name,960,0.0,120.0,1200,0.0,30.0);

    sprintf(name,"xVq2B4[%d]",i);
    xVq2B4[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    sprintf(name,"xVq2Aft[%d]",i);
    xVq2Aft[i] = new TH2F(name,name,200,-3,1,200,-2,2);
    
    sprintf(name,"enuVeshwB4[%d]",i);
    enuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuVeshwAft[%d]",i);
    enuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"emuVeshwB4[%d]",i);
    emuVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"emuVeshwAft[%d]",i);
    emuVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"enuqeVeshwB4[%d]",i);
    enuqeVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);
    sprintf(name,"enuqeVeshwAft[%d]",i);
    enuqeVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,960,0,120);

    sprintf(name,"q2qeVeshwB4[%d]",i);
    q2qeVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2qeVeshwAft[%d]",i);
    q2qeVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);

    sprintf(name,"q2VeshwB4[%d]",i);
    q2VeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);
    sprintf(name,"q2VeshwAft[%d]",i);
    q2VeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,40.0);

    sprintf(name,"xVeshwB4[%d]",i);
    xVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);
    sprintf(name,"xVeshwAft[%d]",i);
    xVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,1200,0.0,30.0);

    sprintf(name,"yVeshwB4[%d]",i);
    yVeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);
    sprintf(name,"yVeshwAft[%d]",i);
    yVeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,200,0.0,1.0);

    sprintf(name,"w2VeshwB4[%d]",i);
    w2VeshwB4[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);
    sprintf(name,"w2VeshwAft[%d]",i);
    w2VeshwAft[i] = new TH2F(name,name,1200,0.0,60.0,300,-5.0,10.0);

  }
  
  ///////////////////////////////////////////////////

  dir->cd();

  Double_t potTGT = 0.0;
  Double_t pot101 = 0.0;

  Int_t nSnarl = 0;
  Bool_t fstSnarl = true;
  Int_t crntSnarl = 0;
  Int_t nRun = 0;
  Int_t crntRun = 0;

  for(Long64_t jentry=0;jentry<nentries;jentry++){
    Long64_t ientry = LoadTree(jentry);
    if(ientry<0) break;
    nb = fChain->GetEntry(jentry);
    nBytes += nb;

    // Beam cuts and POT counting
    
    // Beam configurations: 
    // 0 = unknown, 1 = LE-10, 2 = ME, 3 = HE, 4 = pME, 5 = pHE 

    if(beamcnf!=beamConfig) continue;

    // Exclude special horn current runs

    if(run==8218) continue;
    if(run==8219) continue;
    if(run==8220) continue;
    if(run==8221) continue;
    if(run==8222) continue;
    if(run==8223) continue;
    if(run==8224) continue;
    if(run==8225) continue;
    if(run==8226) continue;
    if(run==8227) continue;
    if(run==8228) continue;
    if(run==8229) continue;
    if(run==8230) continue;
    if(run==8231) continue;
    if(run==8232) continue;
    if(run==8233) continue;
    if(run==8234) continue;
    if(run==8235) continue;
    if(run==8236) continue;
    if(run==8237) continue;
    if(run==8238) continue;

    // Good beam cuts

    if((tor101<0.1||tortgt<0.1||hbw>2.9||vbw>2.9||
	hpos2<-2||hpos2>0||vpos2<0||vpos2>2||
	fabs(closestspill)>2||fabs(hornI)<50)){
      continue;
    }
    
    if(event==0){
      potTGT += tortgt;
      pot101 += tor101;
    }

    // Run and snarl book keeping

    if(fstSnarl){
      crntSnarl = snarl;
      crntRun = run;
      nSnarl++;
      nRun++;
      fstSnarl = false;
    }
    else {
      if(crntRun!=run){
	crntSnarl = snarl;
	crntRun = run;
	nSnarl++;
	nRun++;
      }
      else if(crntSnarl!=snarl){
	crntSnarl = snarl;
	nSnarl++;
      }
    }

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

    // Good event cuts
      
    if(reco_enu>=120) continue;
    // Pitt ND fid vol etc.
    if(pass!=1) continue;
    // Track fit momentum error cut
    if(emu_meth==1 && fabs(trkeqp/trkqp)>0.3) continue;
    // Anti-neutrino cut
    if(trkqp>0) continue;
    // Rock muon cut
    if(reco_trk_vtxz<0.5) continue;
    // Runt event cut
    if((reco_enu - nvtp)/(reco_enu + nvtp)<0.95 && (nvtt*1e9)<50) continue;

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

    // Best approximation of 'standard' CC group good event cuts
      /*
    if(ntrack==0) continue; 
    if(reco_trk_vtxz<=1 || reco_trk_vtxz>=5) continue;
    if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))>=1) continue;
    if(pitt_trk_qual!=1) continue;
    if(reco_dircosneu<=0.6) continue;
    if(trkqp>0) continue;
      */

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

    // Do CC PID cut?
    if(selectCC){
      if(dave_cc_pid<-0.1) continue;
    }

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

    // Fiducial Sample Cuts
      
    // 'Anulus'
    //if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))<0.4 || 
    //   (((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))>=100.0) continue; 
     
    // Quadrants with beam spot 'chopped out'
    //if(reco_trk_vtxx>1.4885) continue;    
    //if(reco_trk_vtxy<=0.0) continue;  
    //if((((reco_trk_vtxx-1.4885)*(reco_trk_vtxx-1.4885))+((reco_trk_vtxy-0.1397)*(reco_trk_vtxy-0.1397)))<0.04) continue; 

    // Z chunks
    if(reco_trk_vtxz<2.5 || reco_trk_vtxz>=3.5) continue;

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

    // Fill histograms before sample selection
    
    enuB4[0]->Fill(reco_enu);
    emuB4[0]->Fill(reco_emu);
    eshwB4[0]->Fill(reco_eshw);
    enuqeB4[0]->Fill(reco_qe_enu);
    thetamuB4[0]->Fill(acos(reco_dircosneu));
    q2qeB4[0]->Fill(reco_qe_q2);
    q2B4[0]->Fill(reco_Q2);
    xB4[0]->Fill(reco_x);
    yB4[0]->Fill(reco_y);
    w2B4[0]->Fill(reco_W2);

    thetamuVemuB4[0]->Fill(reco_emu,acos(reco_dircosneu));
    thetamuVeshwB4[0]->Fill(reco_eshw,acos(reco_dircosneu));
    yVenuB4[0]->Fill(reco_enu,reco_y);
    q2VenuB4[0]->Fill(reco_enu,reco_Q2);
    w2VenuB4[0]->Fill(reco_enu,reco_W2);
    xVenuB4[0]->Fill(reco_enu,reco_x);
    xVq2B4[0]->Fill(log10(reco_x),log10(reco_Q2));
    
    trkVtxB4->Fill(reco_trk_vtxx,reco_trk_vtxy);

    enuVeshwB4[0]->Fill(reco_eshw,reco_enu);
    emuVeshwB4[0]->Fill(reco_eshw,reco_emu);
    enuqeVeshwB4[0]->Fill(reco_eshw,reco_qe_enu);
    q2qeVeshwB4[0]->Fill(reco_eshw,reco_qe_q2);
    q2VeshwB4[0]->Fill(reco_eshw,reco_Q2);
    xVeshwB4[0]->Fill(reco_eshw,reco_x);
    yVeshwB4[0]->Fill(reco_eshw,reco_y);
    w2VeshwB4[0]->Fill(reco_eshw,reco_W2);

    if(reco_Q2>1 && reco_W2>4){
      enuB4[1]->Fill(reco_enu);
      emuB4[1]->Fill(reco_emu);
      eshwB4[1]->Fill(reco_eshw);
      enuqeB4[1]->Fill(reco_qe_enu);
      thetamuB4[1]->Fill(acos(reco_dircosneu));
      q2qeB4[1]->Fill(reco_qe_q2);
      q2B4[1]->Fill(reco_Q2);
      xB4[1]->Fill(reco_x);
      yB4[1]->Fill(reco_y);
      w2B4[1]->Fill(reco_W2);
      
      thetamuVemuB4[1]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwB4[1]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuB4[1]->Fill(reco_enu,reco_y);
      q2VenuB4[1]->Fill(reco_enu,reco_Q2);
      w2VenuB4[1]->Fill(reco_enu,reco_W2);
      xVenuB4[1]->Fill(reco_enu,reco_x);
      xVq2B4[1]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwB4[1]->Fill(reco_eshw,reco_enu);
      emuVeshwB4[1]->Fill(reco_eshw,reco_emu);
      enuqeVeshwB4[1]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwB4[1]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwB4[1]->Fill(reco_eshw,reco_Q2);
      xVeshwB4[1]->Fill(reco_eshw,reco_x);
      yVeshwB4[1]->Fill(reco_eshw,reco_y);
      w2VeshwB4[1]->Fill(reco_eshw,reco_W2);
    }
    if(reco_Q2<=1 && reco_W2>4){
      enuB4[2]->Fill(reco_enu);
      emuB4[2]->Fill(reco_emu);
      eshwB4[2]->Fill(reco_eshw);
      enuqeB4[2]->Fill(reco_qe_enu);
      thetamuB4[2]->Fill(acos(reco_dircosneu));
      q2qeB4[2]->Fill(reco_qe_q2);
      q2B4[2]->Fill(reco_Q2);
      xB4[2]->Fill(reco_x);
      yB4[2]->Fill(reco_y);
      w2B4[2]->Fill(reco_W2);
      
      thetamuVemuB4[2]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwB4[2]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuB4[2]->Fill(reco_enu,reco_y);
      q2VenuB4[2]->Fill(reco_enu,reco_Q2);
      w2VenuB4[2]->Fill(reco_enu,reco_W2);
      xVenuB4[2]->Fill(reco_enu,reco_x);
      xVq2B4[2]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwB4[2]->Fill(reco_eshw,reco_enu);
      emuVeshwB4[2]->Fill(reco_eshw,reco_emu);
      enuqeVeshwB4[2]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwB4[2]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwB4[2]->Fill(reco_eshw,reco_Q2);
      xVeshwB4[2]->Fill(reco_eshw,reco_x);
      yVeshwB4[2]->Fill(reco_eshw,reco_y);
      w2VeshwB4[2]->Fill(reco_eshw,reco_W2);
    }
    if(reco_W2<1.69){
      enuB4[3]->Fill(reco_enu);
      emuB4[3]->Fill(reco_emu);
      eshwB4[3]->Fill(reco_eshw);
      enuqeB4[3]->Fill(reco_qe_enu);
      thetamuB4[3]->Fill(acos(reco_dircosneu));
      q2qeB4[3]->Fill(reco_qe_q2);
      q2B4[3]->Fill(reco_Q2);
      xB4[3]->Fill(reco_x);
      yB4[3]->Fill(reco_y);
      w2B4[3]->Fill(reco_W2);
      
      thetamuVemuB4[3]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwB4[3]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuB4[3]->Fill(reco_enu,reco_y);
      q2VenuB4[3]->Fill(reco_enu,reco_Q2);
      w2VenuB4[3]->Fill(reco_enu,reco_W2);
      xVenuB4[3]->Fill(reco_enu,reco_x);
      xVq2B4[3]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwB4[3]->Fill(reco_eshw,reco_enu);
      emuVeshwB4[3]->Fill(reco_eshw,reco_emu);
      enuqeVeshwB4[3]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwB4[3]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwB4[3]->Fill(reco_eshw,reco_Q2);
      xVeshwB4[3]->Fill(reco_eshw,reco_x);
      yVeshwB4[3]->Fill(reco_eshw,reco_y);
      w2VeshwB4[3]->Fill(reco_eshw,reco_W2);
    }
    if(reco_W2>=1.69 && reco_W2<=4){
      enuB4[4]->Fill(reco_enu);
      emuB4[4]->Fill(reco_emu);
      eshwB4[4]->Fill(reco_eshw);
      enuqeB4[4]->Fill(reco_qe_enu);
      thetamuB4[4]->Fill(acos(reco_dircosneu));
      q2qeB4[4]->Fill(reco_qe_q2);
      q2B4[4]->Fill(reco_Q2);
      xB4[4]->Fill(reco_x);
      yB4[4]->Fill(reco_y);
      w2B4[4]->Fill(reco_W2);
      
      thetamuVemuB4[4]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwB4[4]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuB4[4]->Fill(reco_enu,reco_y);
      q2VenuB4[4]->Fill(reco_enu,reco_Q2);
      w2VenuB4[4]->Fill(reco_enu,reco_W2);
      xVenuB4[4]->Fill(reco_enu,reco_x);
      xVq2B4[4]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwB4[4]->Fill(reco_eshw,reco_enu);
      emuVeshwB4[4]->Fill(reco_eshw,reco_emu);
      enuqeVeshwB4[4]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwB4[4]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwB4[4]->Fill(reco_eshw,reco_Q2);
      xVeshwB4[4]->Fill(reco_eshw,reco_x);
      yVeshwB4[4]->Fill(reco_eshw,reco_y);
      w2VeshwB4[4]->Fill(reco_eshw,reco_W2);
    }

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

    // Select sample

    if(reco_eshw>0.0) continue;

    //if(med_qe_pass!=1) continue;

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

    // Fill histograms after sample selection

    enuAft[0]->Fill(reco_enu);
    emuAft[0]->Fill(reco_emu);
    eshwAft[0]->Fill(reco_eshw);
    enuqeAft[0]->Fill(reco_qe_enu);
    thetamuAft[0]->Fill(acos(reco_dircosneu));
    q2qeAft[0]->Fill(reco_qe_q2);
    q2Aft[0]->Fill(reco_Q2);
    xAft[0]->Fill(reco_x);
    yAft[0]->Fill(reco_y);
    w2Aft[0]->Fill(reco_W2);

    thetamuVemuAft[0]->Fill(reco_emu,acos(reco_dircosneu));
    thetamuVeshwAft[0]->Fill(reco_eshw,acos(reco_dircosneu));
    yVenuAft[0]->Fill(reco_enu,reco_y);
    q2VenuAft[0]->Fill(reco_enu,reco_Q2);
    w2VenuAft[0]->Fill(reco_enu,reco_W2);
    xVenuAft[0]->Fill(reco_enu,reco_x);
    xVq2Aft[0]->Fill(log10(reco_x),log10(reco_Q2));
    
    trkVtxAft->Fill(reco_trk_vtxx,reco_trk_vtxy);

    enuVeshwAft[0]->Fill(reco_eshw,reco_enu);
    emuVeshwAft[0]->Fill(reco_eshw,reco_emu);
    enuqeVeshwAft[0]->Fill(reco_eshw,reco_qe_enu);
    q2qeVeshwAft[0]->Fill(reco_eshw,reco_qe_q2);
    q2VeshwAft[0]->Fill(reco_eshw,reco_Q2);
    xVeshwAft[0]->Fill(reco_eshw,reco_x);
    yVeshwAft[0]->Fill(reco_eshw,reco_y);
    w2VeshwAft[0]->Fill(reco_eshw,reco_W2);

    if(reco_Q2>1 && reco_W2>4){
      enuAft[1]->Fill(reco_enu);
      emuAft[1]->Fill(reco_emu);
      eshwAft[1]->Fill(reco_eshw);
      enuqeAft[1]->Fill(reco_qe_enu);
      thetamuAft[1]->Fill(acos(reco_dircosneu));
      q2qeAft[1]->Fill(reco_qe_q2);
      q2Aft[1]->Fill(reco_Q2);
      xAft[1]->Fill(reco_x);
      yAft[1]->Fill(reco_y);
      w2Aft[1]->Fill(reco_W2);
      
      thetamuVemuAft[1]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwAft[1]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuAft[1]->Fill(reco_enu,reco_y);
      q2VenuAft[1]->Fill(reco_enu,reco_Q2);
      w2VenuAft[1]->Fill(reco_enu,reco_W2);
      xVenuAft[1]->Fill(reco_enu,reco_x);
      xVq2Aft[1]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwAft[1]->Fill(reco_eshw,reco_enu);
      emuVeshwAft[1]->Fill(reco_eshw,reco_emu);
      enuqeVeshwAft[1]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwAft[1]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwAft[1]->Fill(reco_eshw,reco_Q2);
      xVeshwAft[1]->Fill(reco_eshw,reco_x);
      yVeshwAft[1]->Fill(reco_eshw,reco_y);
      w2VeshwAft[1]->Fill(reco_eshw,reco_W2);
    }
    if(reco_Q2<=1 && reco_W2>4){
      enuAft[2]->Fill(reco_enu);
      emuAft[2]->Fill(reco_emu);
      eshwAft[2]->Fill(reco_eshw);
      enuqeAft[2]->Fill(reco_qe_enu);
      thetamuAft[2]->Fill(acos(reco_dircosneu));
      q2qeAft[2]->Fill(reco_qe_q2);
      q2Aft[2]->Fill(reco_Q2);
      xAft[2]->Fill(reco_x);
      yAft[2]->Fill(reco_y);
      w2Aft[2]->Fill(reco_W2);
      
      thetamuVemuAft[2]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwAft[2]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuAft[2]->Fill(reco_enu,reco_y);
      q2VenuAft[2]->Fill(reco_enu,reco_Q2);
      w2VenuAft[2]->Fill(reco_enu,reco_W2);
      xVenuAft[2]->Fill(reco_enu,reco_x);
      xVq2Aft[2]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwAft[2]->Fill(reco_eshw,reco_enu);
      emuVeshwAft[2]->Fill(reco_eshw,reco_emu);
      enuqeVeshwAft[2]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwAft[2]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwAft[2]->Fill(reco_eshw,reco_Q2);
      xVeshwAft[2]->Fill(reco_eshw,reco_x);
      yVeshwAft[2]->Fill(reco_eshw,reco_y);
      w2VeshwAft[2]->Fill(reco_eshw,reco_W2);
    }
    if(reco_W2<1.69){
      enuAft[3]->Fill(reco_enu);
      emuAft[3]->Fill(reco_emu);
      eshwAft[3]->Fill(reco_eshw);
      enuqeAft[3]->Fill(reco_qe_enu);
      thetamuAft[3]->Fill(acos(reco_dircosneu));
      q2qeAft[3]->Fill(reco_qe_q2);
      q2Aft[3]->Fill(reco_Q2);
      xAft[3]->Fill(reco_x);
      yAft[3]->Fill(reco_y);
      w2Aft[3]->Fill(reco_W2);
      
      thetamuVemuAft[3]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwAft[3]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuAft[3]->Fill(reco_enu,reco_y);
      q2VenuAft[3]->Fill(reco_enu,reco_Q2);
      w2VenuAft[3]->Fill(reco_enu,reco_W2);
      xVenuAft[3]->Fill(reco_enu,reco_x);
      xVq2Aft[3]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwAft[3]->Fill(reco_eshw,reco_enu);
      emuVeshwAft[3]->Fill(reco_eshw,reco_emu);
      enuqeVeshwAft[3]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwAft[3]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwAft[3]->Fill(reco_eshw,reco_Q2);
      xVeshwAft[3]->Fill(reco_eshw,reco_x);
      yVeshwAft[3]->Fill(reco_eshw,reco_y);
      w2VeshwAft[3]->Fill(reco_eshw,reco_W2);
    }
    if(reco_W2>=1.69 && reco_W2<=4){
      enuAft[4]->Fill(reco_enu);
      emuAft[4]->Fill(reco_emu);
      eshwAft[4]->Fill(reco_eshw);
      enuqeAft[4]->Fill(reco_qe_enu);
      thetamuAft[4]->Fill(acos(reco_dircosneu));
      q2qeAft[4]->Fill(reco_qe_q2);
      q2Aft[4]->Fill(reco_Q2);
      xAft[4]->Fill(reco_x);
      yAft[4]->Fill(reco_y);
      w2Aft[4]->Fill(reco_W2);
      
      thetamuVemuAft[4]->Fill(reco_emu,acos(reco_dircosneu));
      thetamuVeshwAft[4]->Fill(reco_eshw,acos(reco_dircosneu));
      yVenuAft[4]->Fill(reco_enu,reco_y);
      q2VenuAft[4]->Fill(reco_enu,reco_Q2);
      w2VenuAft[4]->Fill(reco_enu,reco_W2);
      xVenuAft[4]->Fill(reco_enu,reco_x);
      xVq2Aft[4]->Fill(log10(reco_x),log10(reco_Q2));

      enuVeshwAft[4]->Fill(reco_eshw,reco_enu);
      emuVeshwAft[4]->Fill(reco_eshw,reco_emu);
      enuqeVeshwAft[4]->Fill(reco_eshw,reco_qe_enu);
      q2qeVeshwAft[4]->Fill(reco_eshw,reco_qe_q2);
      q2VeshwAft[4]->Fill(reco_eshw,reco_Q2);
      xVeshwAft[4]->Fill(reco_eshw,reco_x);
      yVeshwAft[4]->Fill(reco_eshw,reco_y);
      w2VeshwAft[4]->Fill(reco_eshw,reco_W2);
    }

  }

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

  // Output histograms to file

  for(Int_t i=0;i<5;i++){
      
    enuB4[i]->SetDirectory(fout);
    enuAft[i]->SetDirectory(fout);
    
    emuB4[i]->SetDirectory(fout);
    emuAft[i]->SetDirectory(fout);
      
    eshwB4[i]->SetDirectory(fout);
    eshwAft[i]->SetDirectory(fout);
      
    enuqeB4[i]->SetDirectory(fout);
    enuqeAft[i]->SetDirectory(fout);
      
    thetamuB4[i]->SetDirectory(fout);
    thetamuAft[i]->SetDirectory(fout);
      
    q2qeB4[i]->SetDirectory(fout);
    q2qeAft[i]->SetDirectory(fout);
      
    q2B4[i]->SetDirectory(fout);
    q2Aft[i]->SetDirectory(fout);
      
    xB4[i]->SetDirectory(fout);
    xAft[i]->SetDirectory(fout);
      
    yB4[i]->SetDirectory(fout);
    yAft[i]->SetDirectory(fout);
      
    w2B4[i]->SetDirectory(fout);
    w2Aft[i]->SetDirectory(fout);

    thetamuVemuB4[i]->SetDirectory(fout);
    thetamuVemuAft[i]->SetDirectory(fout);
      
    thetamuVeshwB4[i]->SetDirectory(fout);
    thetamuVeshwAft[i]->SetDirectory(fout);
      
    yVenuB4[i]->SetDirectory(fout);
    yVenuAft[i]->SetDirectory(fout);
      
    q2VenuB4[i]->SetDirectory(fout);
    q2VenuAft[i]->SetDirectory(fout);

    w2VenuB4[i]->SetDirectory(fout);
    w2VenuAft[i]->SetDirectory(fout);

    xVenuB4[i]->SetDirectory(fout);
    xVenuAft[i]->SetDirectory(fout);
      
    xVq2B4[i]->SetDirectory(fout);
    xVq2Aft[i]->SetDirectory(fout);
    
    enuVeshwB4[i]->SetDirectory(fout);
    enuVeshwAft[i]->SetDirectory(fout);
    
    emuVeshwB4[i]->SetDirectory(fout);
    emuVeshwAft[i]->SetDirectory(fout);
      
    enuqeVeshwB4[i]->SetDirectory(fout);
    enuqeVeshwAft[i]->SetDirectory(fout);
      
    q2qeVeshwB4[i]->SetDirectory(fout);
    q2qeVeshwAft[i]->SetDirectory(fout);
      
    q2VeshwB4[i]->SetDirectory(fout);
    q2VeshwAft[i]->SetDirectory(fout);
      
    xVeshwB4[i]->SetDirectory(fout);
    xVeshwAft[i]->SetDirectory(fout);
      
    yVeshwB4[i]->SetDirectory(fout);
    yVeshwAft[i]->SetDirectory(fout);
      
    w2VeshwB4[i]->SetDirectory(fout);
    w2VeshwAft[i]->SetDirectory(fout);
  
  }
  
  trkVtxB4->SetDirectory(fout);
  trkVtxAft->SetDirectory(fout);

  fout->Write();
 
  cout << "Number of snarls:   " << nSnarl << endl;
  cout << "Number of runs:     " << nRun << endl;

  cout << "POTs from TorTGT:   " << potTGT << "*e12" << endl; 
  cout << "POTs from Tor101:   " << pot101 << "*e12" << endl; 

}

//***************************************************
//***************************************************
//***************************************************