#ifndef newmadqeid_cxx
#define newmadqeid_cxx

#include "Mad/NewMadQEID.h"
#include "Mad/NearbyEvents.h"
#include "Mad/HoughTrans.h"
#include "Mad/MadBase.h"
#include <cmath>
#include "TFile.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TVector3.h"
#include "TF1.h"

#define PI 3.14159265

//**********************************************************
NewMadQEID::NewMadQEID() 
  : nBins(0), maxE(20.0),
    nShows(-1), useNshows(true),
    houghPeak(-1), useHoughPeak(true),
    recoWsqd(-100.0), useRecoWsqd(true),
    rePHfrac(-1.0), useRePHfrac(true),
    nVHShits(-1), useNVHShits(true),
    nSubShows(-1), useNsubshows(false),
    vhsPH(-1.0), useVhsPH(false),
    qeKinEnuDif(-1.0), useQeKinEnuDif(false),
    userVar(-1.0), useUserVariable(false),
    pdfsNormed(false),
    binningSet(false),
    protonAngU(-10.0),
    protonAngV(-10.0),
    stdHepProtAngU(-10.0),
    stdHepProtAngV(-10.0),
    allHitAngU(-10.0),
    allHitAngV(-10.0),
    VHSangU(-10.0),
    VHSangV(-10.0),
    protonUcomp(-1.0),
    protonVcomp(-1.0),
    protonZcomp(-1.0),
    protonEfromEvent(-1.0),
    protonEfromMassMom(-1.0),
    stdHepProtUcomp(-1.0),
    stdHepProtVcomp(-1.0),
    stdHepProtZcomp(-1.0),
    stdHepProtE(-1.0),
    muonUcomp(-1.0),
    muonVcomp(-1.0),
    muonZcomp(-1.0),
    muonE(-1.0),
    stdHepMuUcomp(-1.0),
    stdHepMuVcomp(-1.0),
    stdHepMuZcomp(-1.0),
    stdHepMuE(-1.0),
    houghRmsU(-1.0),
    houghRmsV(-1.0),
    phIn0_5sigCone(-1.0),
    phIn1_0sigCone(-1.0),
    phIn1_5sigCone(-1.0),
    phIn2_0sigCone(-1.0),
    teIn0_5sigCone(-1.0),
    teIn1_0sigCone(-1.0),
    teIn1_5sigCone(-1.0),
    teIn2_0sigCone(-1.0),
    teTot(-1.0),
    nhIn0_5sigCone(-1),
    nhIn1_0sigCone(-1),
    nhIn1_5sigCone(-1),
    nhIn2_0sigCone(-1),
    rePH(-1.0),
    totPH(-1.0),
    phInProtR(-1.0),
    shwAngSpread(-10.0)

{

  for(Int_t i=0;i<30;i++){

    nShowQE[i] = 0;
    nShowBG[i] = 0;
    hPeakQE[i] = 0;
    hPeakBG[i] = 0;
    recoW2QE[i] = 0;
    recoW2BG[i] = 0;
    rePHfracQE[i] = 0;
    rePHfracBG[i] = 0;
    nVHShitsQE[i] = 0;
    nVHShitsBG[i] = 0;
    nSSQE[i] = 0;
    nSSBG[i] = 0;
    vhsPHQE[i] = 0;
    vhsPHBG[i] = 0;
    qeKinEnuDifQE[i] = 0;
    qeKinEnuDifBG[i] = 0;
    userVarQE[i] = 0;
    userVarBG[i] = 0;
    trainingHists[i] = 0;

  }

  cutsHist = 0;
  nBinsHist = 0;
  eBins = 0;

  InitPrange();

}

//********************************************************
void NewMadQEID::SetEnergyBinning(Float_t axArray[31])
{

  if(binningSet) return;

  eBins = new TH1F("eBins","eBins",960,0,120);

  //cout << "Setting energy binning..." << endl;

  Int_t count = 0;
  Float_t maxVal = 0.0;
  for(Int_t i=0;i<31;i++){
    if(i==0){
      binArray[i] = axArray[i];
      eBins->Fill(axArray[i]+0.01); // + 0.01 keeps the correct binning when reading in PDFs
      //cout << "Bin " << i << " has low edge " << binArray[i]  << endl;
    }
    else if(i>0 && (axArray[i]>axArray[i-1])){
      count++;
      binArray[i] = axArray[i];
      maxVal = binArray[i];
      eBins->Fill(axArray[i]+0.01);
      //cout << "Bin " << i << " has low edge " << binArray[i]  << endl;
    }
    else {
      binArray[i] = 0.0;
      //cout << "Bin " << i << " has low edge " << binArray[i]  << endl;
    }
  }

  nBins = count;

  //cout << "nBins: " << nBins << endl;

  maxE = maxVal;

  //cout << "maxE:  " << maxE << endl;

  nBinsHist = new TH1F("nBinsHist","nBinsHist",30,0,30);
  nBinsHist->Fill(nBins);
  binningSet = true;

  return;

}

//********************************************************
void NewMadQEID::InitPrange()
{

  prange1 = new TF1("prange1","0.1083-1.3817*x-38.5143*x*x+279.1245*x*x*x-124.3302*x*x*x*x+17.2100*x*x*x*x*x",0.0,2.7);
  prange2 = new TF1("prange2","-777.5115+691.6554*x-3.1098*x*x",2.7,11);

}

//********************************************************
void NewMadQEID::InitPDFs()
{

  TH1::AddDirectory(false); // do not assign histograms to TDirectory

  for(Int_t i=0;i<30;i++){
    char name[30];

    if(useRePHfrac){
      sprintf(name,"rePHfracQE[%d]",i);
      rePHfracQE[i] = new TH1F(name,name,25,0.0,1.0);
      sprintf(name,"rePHfracBG[%d]",i);
      rePHfracBG[i] = new TH1F(name,name,25,0.0,1.0);
    }

    if(useRecoWsqd){
      sprintf(name,"recoW2QE[%d]",i);
      recoW2QE[i] = new TH1F(name,name,75,-5,10);
      sprintf(name,"recoW2BG[%d]",i);
      recoW2BG[i] = new TH1F(name,name,75,-5,10);
    }

    if(useNshows){
      sprintf(name,"nShowQE[%d]",i);
      nShowQE[i] = new TH1F(name,name,10,0,10);
      sprintf(name,"nShowBG[%d]",i);
      nShowBG[i] = new TH1F(name,name,10,0,10);
    }

    if(useHoughPeak){
      sprintf(name,"hPeakQE[%d]",i);
      hPeakQE[i] = new TH1F(name,name,50,0,50);
      sprintf(name,"hPeakBG[%d]",i);
      hPeakBG[i] = new TH1F(name,name,50,0,50);
    }

    if(useNVHShits){
      sprintf(name,"nVHShitsQE[%d]",i);
      nVHShitsQE[i] = new TH1F(name,name,50,0,50);
      sprintf(name,"nVHShitsBG[%d]",i);
      nVHShitsBG[i] = new TH1F(name,name,50,0,50);
    }

    if(useNsubshows){
      sprintf(name,"nSSQE[%d]",i);
      nSSQE[i] = new TH1F(name,name,50,0,50);
      sprintf(name,"nSSBG[%d]",i);
      nSSBG[i] = new TH1F(name,name,50,0,50);
    }

    if(useVhsPH){
      sprintf(name,"vhsPHQE[%d]",i);
      vhsPHQE[i] = new TH1F(name,name,50,0,500000);
      sprintf(name,"vhsPHBG[%d]",i);
      vhsPHBG[i] = new TH1F(name,name,50,0,500000);
    }

    if(useQeKinEnuDif){
      sprintf(name,"qeKinEnuDifQE[%d]",i);
      qeKinEnuDifQE[i] = new TH1F(name,name,480,-120,120);
      sprintf(name,"qeKinEnuDifBG[%d]",i);
      qeKinEnuDifBG[i] = new TH1F(name,name,480,-120,120);
    }

    if(useUserVariable){
      sprintf(name,"userVarQE[%d]",i);
      userVarQE[i] = new TH1F(name,name,50,0,1);
      sprintf(name,"userVarBG[%d]",i);
      userVarBG[i] = new TH1F(name,name,50,0,1);
    }

  }  

  TH1::AddDirectory(kTRUE); // revert above

}

//********************************************************
void NewMadQEID::InitTrainingHists()
{

  TH1::AddDirectory(false); // do not assign histograms to TDirectory

  for(Int_t i=0;i<30;i++){
    char name[30];
    sprintf(name,"trainingHists[%d]",i);
    trainingHists[i] = new TH1F(name,name,1200,-6.0,6.0);
  }

  Float_t axArray[31];

  if(!binningSet){
    eBins = new TH1F("eBins","eBins",960,0,120);
    Float_t tmpAxArray[19] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,
			      5.0,6.0,7.0,8.0,10.0,12.0,14.0,17.0,20.0};
    for(Int_t i=0;i<19;i++){
      binArray[i] = tmpAxArray[i];
      eBins->Fill(tmpAxArray[i]);
    }
    for(Int_t i=19;i<31;i++){
      binArray[i] = 0.0;
    }
    for(Int_t i=0;i<19;i++){
      axArray[i] = binArray[i];
    }
    for(Int_t i=19;i<31;i++){
      axArray[i] = (axArray[18]+((i+1)*0.1));
    }
    nBins = 18;
    maxE = 20.0;
    nBinsHist = new TH1F("nBinsHist","nBinsHist",30,0,30);
    nBinsHist->Fill(nBins);
    binningSet = true;
  }

  else {
    //Float_t tmpaxArray[nBins];
    for(Int_t i=0;i<(nBins+1);i++){
      axArray[i] = binArray[i];
    }
    for(Int_t i=(nBins+1);i<31;i++){
      axArray[i] = (binArray[nBins]+((i+1)*0.1));
    }
  }

  //cout << "Creating cuts histo." << endl;
  //for(Int_t i=0;i<31;i++){
  //  cout << "Cuts histo array value " << i << " is " << axArray[i]  << endl;
  //}
  // cutsHist = new TH1F("cuthist","cuthist",nBins,axArray);
  cutsHist = new TH1F("cuthist","cuthist",30,axArray);
  //cout << "Cuts hists created." << endl;

  TH1::AddDirectory(true); // revert above

  return;

}

//********************************************************
void NewMadQEID::ResetPDFs()
{

  for(Int_t i=0;i<30;i++){

    if(useRePHfrac){
      rePHfracQE[i]->Reset();
      rePHfracBG[i]->Reset();
    }
    if(useRecoWsqd){
      recoW2QE[i]->Reset();
      recoW2BG[i]->Reset();
    }
    if(useNshows){
      nShowQE[i]->Reset();
      nShowBG[i]->Reset();
    }
    if(useHoughPeak){
      hPeakQE[i]->Reset();
      hPeakBG[i]->Reset();
    }
    if(useNVHShits){
      nVHShitsQE[i]->Reset();
      nVHShitsBG[i]->Reset();
    }
    if(useNsubshows){
      nSSQE[i]->Reset();
      nSSBG[i]->Reset();
    }
    if(useVhsPH){
      vhsPHQE[i]->Reset();
      vhsPHBG[i]->Reset();
    }
    if(useQeKinEnuDif){
      qeKinEnuDifQE[i]->Reset();
      qeKinEnuDifBG[i]->Reset();
    }
    if(useUserVariable){
      userVarQE[i]->Reset();
      userVarBG[i]->Reset();
    }
  }

}

//********************************************************
void NewMadQEID::ReInitVars()
{

  nShows = -1;
  houghPeak = -1;
  recoWsqd = -100.0;
  rePHfrac = -1.0;
  nVHShits = -1;
  nSubShows = -1;
  vhsPH = -1.0;
  qeKinEnuDif = -1.0;
  userVar = -1.0; 
  
  protonAngU = -10.0;
  protonAngV = -10.0;
  stdHepProtAngU = -10.0;
  stdHepProtAngV = -10.0;
  allHitAngU = -10.0;
  allHitAngV = -10.0;
  VHSangU = -10.0;
  VHSangV = -10.0;
  protonUcomp = -1.0;
  protonVcomp = -1.0;
  protonZcomp = -1.0;
  protonEfromEvent = -1.0;
  protonEfromMassMom = -1.0;
  stdHepProtUcomp = -1.0;
  stdHepProtVcomp = -1.0;
  stdHepProtZcomp = -1.0;
  stdHepProtE = -1.0;

  muonUcomp = -1.0;
  muonVcomp = -1.0;
  muonZcomp = -1.0;
  muonE = -1.0;
  stdHepMuUcomp = -1.0;
  stdHepMuVcomp = -1.0;
  stdHepMuZcomp = -1.0;
  stdHepMuE = -1.0;

  houghRmsU = -1.0;
  houghRmsV = -1.0;

  phIn0_5sigCone = -1.0;
  phIn1_0sigCone = -1.0;
  phIn1_5sigCone = -1.0;
  phIn2_0sigCone = -1.0;

  teIn0_5sigCone = -1.0;
  teIn1_0sigCone = -1.0;
  teIn1_5sigCone = -1.0;
  teIn2_0sigCone = -1.0;

  teTot = -1.0;

  nhIn0_5sigCone = -1;
  nhIn1_0sigCone = -1;
  nhIn1_5sigCone = -1;
  nhIn2_0sigCone = -1;

  rePH = -1.0;
  totPH = -1.0;

  phInProtR = -1.0;
  shwAngSpread = -10.0;

  return;

}

//********************************************************
Bool_t NewMadQEID::CalcQEVars(MadBase* mb, Int_t evtindex,
			      Float_t reco_enu, Float_t reco_emu)
{

  if(!binningSet){
    Float_t axArray[31] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,
			   5.0,6.0,7.0,8.0,10.0,12.0,14.0,17.0,20.0,
			   0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    SetEnergyBinning(axArray);
  }

  ReInitVars();

  if(!mb) return false;
  const NtpSREvent* evt = mb->GetEvent(evtindex);
  if(!evt) return false;

  nShows = evt->nshower;

  Int_t trackindex = -1;
  const NtpSRTrack* trk = mb->GetLargestTrackFromEvent(evtindex,trackindex);
  if(trackindex==-1) return false;

  Int_t trackVtxPlane = trk->vtx.plane;

  // QE kinematics

  TVector3* prot3mom = new TVector3(-100.0,-100.0,-100.0);

  Float_t muonDcosZ = trk->vtx.dcosz;
  Float_t muonDcosU = trk->vtx.dcosu;
  Float_t muonDcosV = trk->vtx.dcosv;
  Float_t muonEnergy = reco_emu;
  Float_t muonMass = 0.10555;
  Float_t muonMom = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
  Float_t nuEnergy = reco_enu;
  Float_t nuDcosZ = cos(-0.058); // assume nu has no x-comp
  Float_t nuDcosY = sin(-0.058); // 3deg in radians
  Float_t nuDcosUV = nuDcosY*(1.0/sqrt(2.0));

  // set 'muon related' data members
  muonUcomp = (muonDcosU*muonMom);
  muonVcomp = (muonDcosV*muonMom);
  muonZcomp = (muonDcosZ*muonMom);
  muonE = reco_emu;

  Float_t dProtMomZ = nuEnergy*nuDcosZ - muonMom*muonDcosZ;
  Float_t dProtMomU = nuEnergy*nuDcosUV - muonMom*muonDcosU;
  Float_t dProtMomV = nuEnergy*nuDcosUV - muonMom*muonDcosV;
  Float_t protMom = sqrt(dProtMomZ*dProtMomZ + dProtMomU*dProtMomU 
			 + dProtMomV*dProtMomV);
  if(protMom>0){
    prot3mom->SetX(dProtMomU);
    prot3mom->SetY(dProtMomV);
    prot3mom->SetZ(dProtMomZ);
  }
  Float_t neutronMass = 0.93957;
  Float_t protMass = 0.93827;
  Float_t protEnergyFromEvent = reco_enu + neutronMass - reco_emu;
  Float_t protEnergyFromMassMom = sqrt(protMom*protMom + protMass*protMass);

  // could use GetDirectionCase() and GetAngle() code to determine
  // the proton direction in the U and V planes but also should test
  // just using the trig functions directly

  Float_t protAngU = GetAngle(prot3mom->X(),prot3mom->Z());
  Float_t protAngV = GetAngle(prot3mom->Y(),prot3mom->Z());

  // set 'proton related' data members
  protonAngU = protAngU;
  protonAngV = protAngV;
  protonUcomp = prot3mom->X();
  protonVcomp = prot3mom->Y();
  protonZcomp = prot3mom->Z();
  protonEfromEvent = protEnergyFromEvent;
  protonEfromMassMom = protEnergyFromMassMom;

  Float_t predProtRange = -1.0;

  if(protonEfromMassMom<2.7 && (prange1->Eval(protEnergyFromEvent-protMass)>=0)){
    predProtRange = prange1->Eval(protEnergyFromEvent-protMass);
  }
  if(protonEfromMassMom>=2.7 && (prange2->Eval(protEnergyFromEvent-protMass)>=0)){
    predProtRange = prange2->Eval(protEnergyFromEvent-protMass);
  }
  if(((prange1->Eval(protEnergyFromEvent-protMass)<0) && protonEfromMassMom<2.7) || 
     ((prange2->Eval(protEnergyFromEvent-protMass)<0) && protonEfromMassMom>=2.7)){
    predProtRange = 0.0;
  }

  Float_t onePlaneIngcm2 = 21.0;
  if(predProtRange!=-1) predProtRange = (predProtRange/onePlaneIngcm2);

  userVar = predProtRange;

  Float_t halfSig = (PI/6.0);

  Float_t ev0_5sigLowU  = protAngU-halfSig;
  Float_t ev0_5sigHighU = protAngU+halfSig;
  Float_t ev1_0sigLowU  = protAngU-(2.0*halfSig);
  Float_t ev1_0sigHighU = protAngU+(2.0*halfSig);
  Float_t ev1_5sigLowU  = protAngU-(3.0*halfSig);
  Float_t ev1_5sigHighU = protAngU+(3.0*halfSig);
  Float_t ev2_0sigLowU  = protAngU-(4.0*halfSig);
  Float_t ev2_0sigHighU = protAngU+(4.0*halfSig);

  Float_t ev0_5sigLowV  = protAngV-halfSig;
  Float_t ev0_5sigHighV = protAngV+halfSig;
  Float_t ev1_0sigLowV  = protAngV-(2.0*halfSig);
  Float_t ev1_0sigHighV = protAngV+(2.0*halfSig);
  Float_t ev1_5sigLowV  = protAngV-(3.0*halfSig);
  Float_t ev1_5sigHighV = protAngV+(3.0*halfSig);
  Float_t ev2_0sigLowV  = protAngV-(4.0*halfSig);
  Float_t ev2_0sigHighV = protAngV+(4.0*halfSig);

  Float_t ph0_5sigU = 0.0;
  Float_t ph1_0sigU = 0.0;
  Float_t ph1_5sigU = 0.0;
  Float_t ph2_0sigU = 0.0;

  Float_t ph0_5sigV = 0.0;
  Float_t ph1_0sigV = 0.0;
  Float_t ph1_5sigV = 0.0;
  Float_t ph2_0sigV = 0.0;

  Int_t nh0_5sigU = 0;
  Int_t nh1_0sigU = 0;
  Int_t nh1_5sigU = 0;
  Int_t nh2_0sigU = 0;

  Int_t nh0_5sigV = 0;
  Int_t nh1_0sigV = 0;
  Int_t nh1_5sigV = 0;
  Int_t nh2_0sigV = 0;

  Float_t te0_5sigU = 0.0;
  Float_t te1_0sigU = 0.0;
  Float_t te1_5sigU = 0.0;
  Float_t te2_0sigU = 0.0;

  Float_t te0_5sigV = 0.0;
  Float_t te1_0sigV = 0.0;
  Float_t te1_5sigV = 0.0;
  Float_t te2_0sigV = 0.0;

  Float_t totStdHepE = 0.0;

  //cout << "----------------" << endl;

  // get and set StdHep 'proton and muon related' data members
  Int_t mc_index = -1;
  const NtpMCTruth* mcTruth = mb->GetTruthForRecoTH(evtindex,mc_index);
  if(mc_index!=-1){ 
    for(Int_t ish=mcTruth->stdhep[0];ish<=mcTruth->stdhep[1];ish++){
      const NtpMCStdHep* stdHep = mb->GetStdHep(ish);
      // proton
      if((stdHep->IdHEP)==2212 && (stdHep->IstHEP)==1){
	//cout << "Got proton." << endl;
	Float_t stdHepPmomX = stdHep->p4[0];
	Float_t stdHepPmomY = stdHep->p4[1];
	Float_t invsqrt2 = (1.0/sqrt(2.0));
	stdHepProtUcomp = ((stdHepPmomX+stdHepPmomY)*invsqrt2);
	stdHepProtVcomp = ((-stdHepPmomX+stdHepPmomY)*invsqrt2);
	stdHepProtZcomp = stdHep->p4[2];
	stdHepProtE = stdHep->p4[3];
	stdHepProtAngU = GetAngle(stdHepProtUcomp,stdHepProtZcomp);
	stdHepProtAngV = GetAngle(stdHepProtVcomp,stdHepProtZcomp);

	totStdHepE+=stdHep->p4[3];

	if(stdHepProtAngV>ev0_5sigLowV && stdHepProtAngV<ev0_5sigHighV){
	  te0_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngV>ev1_0sigLowV && stdHepProtAngV<ev1_0sigHighV){
	  te1_0sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngV>ev1_5sigLowV && stdHepProtAngV<ev1_5sigHighV){
	  te1_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngV>ev2_0sigLowV && stdHepProtAngV<ev2_0sigHighV){
	  te2_0sigV+=(0.5*stdHep->p4[3]);
	}

	if(stdHepProtAngU>ev0_5sigLowU && stdHepProtAngU<ev0_5sigHighU){
	  te0_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngU>ev1_0sigLowU && stdHepProtAngU<ev1_0sigHighU){
	  te1_0sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngU>ev1_5sigLowU && stdHepProtAngU<ev1_5sigHighU){
	  te1_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepProtAngU>ev2_0sigLowU && stdHepProtAngU<ev2_0sigHighU){
	  te2_0sigU+=(0.5*stdHep->p4[3]);
	}

      }
      //muon
      else if((stdHep->IdHEP)==13 && (stdHep->IstHEP)==1){
	Float_t stdHepMuMomX = stdHep->p4[0];
	Float_t stdHepMuMomY = stdHep->p4[1];
	Float_t invsqrt2 = (1.0/sqrt(2.0));
	stdHepMuUcomp = ((stdHepMuMomX+stdHepMuMomY)*invsqrt2);
	stdHepMuVcomp = ((-stdHepMuMomX+stdHepMuMomY)*invsqrt2);
	stdHepMuZcomp = stdHep->p4[2];
	stdHepMuE = stdHep->p4[3];
      }
      //charged pion
      else if((stdHep->IdHEP)==211 && (stdHep->IstHEP)==1){

	//cout << "Got charged pion." << endl;

	Float_t stdHepPimomX = stdHep->p4[0];
	Float_t stdHepPimomY = stdHep->p4[1];
	Float_t invsqrt2 = (1.0/sqrt(2.0));
	Float_t stdHepCpiUcomp = ((stdHepPimomX+stdHepPimomY)*invsqrt2);
	Float_t stdHepCpiVcomp = ((-stdHepPimomX+stdHepPimomY)*invsqrt2);
	Float_t stdHepCpiZcomp = stdHep->p4[2];
	Float_t stdHepCpiAngU = GetAngle(stdHepCpiUcomp,stdHepCpiZcomp);
	Float_t stdHepCpiAngV = GetAngle(stdHepCpiVcomp,stdHepCpiZcomp);

	totStdHepE+=stdHep->p4[3];

	if(stdHepCpiAngV>ev0_5sigLowV && stdHepCpiAngV<ev0_5sigHighV){
	  te0_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngV>ev1_0sigLowV && stdHepCpiAngV<ev1_0sigHighV){
	  te1_0sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngV>ev1_5sigLowV && stdHepCpiAngV<ev1_5sigHighV){
	  te1_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngV>ev2_0sigLowV && stdHepCpiAngV<ev2_0sigHighV){
	  te2_0sigV+=(0.5*stdHep->p4[3]);
	}

	if(stdHepCpiAngU>ev0_5sigLowU && stdHepCpiAngU<ev0_5sigHighU){
	  te0_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngU>ev1_0sigLowU && stdHepCpiAngU<ev1_0sigHighU){
	  te1_0sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngU>ev1_5sigLowU && stdHepCpiAngU<ev1_5sigHighU){
	  te1_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepCpiAngU>ev2_0sigLowU && stdHepCpiAngU<ev2_0sigHighU){
	  te2_0sigU+=(0.5*stdHep->p4[3]);
	}

      }
      //neutral pion
      else if((stdHep->IdHEP)==111 && (stdHep->IstHEP)==1){

	//cout << "Got neutral pion." << endl;

	Float_t stdHepPimomX = stdHep->p4[0];
	Float_t stdHepPimomY = stdHep->p4[1];
	Float_t invsqrt2 = (1.0/sqrt(2.0));
	Float_t stdHepNpiUcomp = ((stdHepPimomX+stdHepPimomY)*invsqrt2);
	Float_t stdHepNpiVcomp = ((-stdHepPimomX+stdHepPimomY)*invsqrt2);
	Float_t stdHepNpiZcomp = stdHep->p4[2];
	Float_t stdHepNpiAngU = GetAngle(stdHepNpiUcomp,stdHepNpiZcomp);
	Float_t stdHepNpiAngV = GetAngle(stdHepNpiVcomp,stdHepNpiZcomp);

	totStdHepE+=stdHep->p4[3];

	if(stdHepNpiAngV>ev0_5sigLowV && stdHepNpiAngV<ev0_5sigHighV){
	  te0_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngV>ev1_0sigLowV && stdHepNpiAngV<ev1_0sigHighV){
	  te1_0sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngV>ev1_5sigLowV && stdHepNpiAngV<ev1_5sigHighV){
	  te1_5sigV+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngV>ev2_0sigLowV && stdHepNpiAngV<ev2_0sigHighV){
	  te2_0sigV+=(0.5*stdHep->p4[3]);
	}

	if(stdHepNpiAngU>ev0_5sigLowU && stdHepNpiAngU<ev0_5sigHighU){
	  te0_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngU>ev1_0sigLowU && stdHepNpiAngU<ev1_0sigHighU){
	  te1_0sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngU>ev1_5sigLowU && stdHepNpiAngU<ev1_5sigHighU){
	  te1_5sigU+=(0.5*stdHep->p4[3]);
	}
	else if(stdHepNpiAngU>ev2_0sigLowU && stdHepNpiAngU<ev2_0sigHighU){
	  te2_0sigU+=(0.5*stdHep->p4[3]);
	}

      }
    }
  }

  teIn0_5sigCone = (te0_5sigU+te0_5sigV);
  teIn1_0sigCone = (te1_0sigU+te1_0sigV);
  teIn1_5sigCone = (te1_5sigU+te1_5sigV);
  teIn2_0sigCone = (te2_0sigU+te2_0sigV);

  teTot = totStdHepE;

  //cout << "Total true E: " << teTot  << endl;
  //cout << " In 0.5 cone: " << teIn0_5sigCone << endl;
  //cout << " In 1.0 cone: " << teIn1_0sigCone << endl;
  //cout << " In 1.5 cone: " << teIn1_5sigCone << endl;
  //cout << " In 2.0 cone: " << teIn2_0sigCone << endl;

  // should I create a radial Hough transform, i.e. where the steps in 
  // gradient give equal sectors of a circle (step changes will be non-linear)?

  static HoughTrans HTUpln(80,-4.0,4.0);
  static HoughTrans HTVpln(80,-4.0,4.0);
  HTUpln.ResetHough();
  HTVpln.ResetHough();
  HTUpln.SetVtxz(trk->vtx.z);
  HTVpln.SetVtxz(trk->vtx.z);

  Float_t totalEventPH = 0.0;
  Float_t remainingPH = 0.0;
  Int_t numVHShits = 0;
  Float_t oneMIPinScor = 545.0;
  //Float_t oneMIPinPEs = 5.45;

  // variables to get PH weighted average angle of VHS hits in U and V
  Float_t allHitSummedAnglesU = 0.0;
  Float_t allHitSummedAnglesV = 0.0;
  Int_t allHitUangleCount = 0;
  Int_t allHitVangleCount = 0;
  Float_t VHShitSummedAnglesU = 0.0;
  Float_t VHShitSummedAnglesV = 0.0;
  Int_t VHShitUangleCount = 0;
  Int_t VHShitVangleCount = 0;

  Int_t numAnglesU = 0;
  Float_t lowAngleU = 0.0;
  Float_t highAngleU = 0.0;
  Int_t numAnglesV = 0;
  Float_t lowAngleV = 0.0;
  Float_t highAngleV = 0.0;

  Float_t phInProtRange = 0.0;

  for(Int_t istp=0;istp<evt->nstrip;istp++){

    int stripindex = evt->stp[istp];
    const NtpSRStrip* stp = mb->GetStrip(stripindex);
    if(!stp) continue;

    Float_t zPos = stp->z;
    Float_t tPos = stp->tpos;
    Float_t hitScor = (stp->ph0.sigcor)+(stp->ph1.sigcor);
    Int_t pln = stp->plane; 

    totalEventPH+=hitScor;

    if(IsVPln(pln)){
      Float_t tmpVcomp = (tPos - trk->vtx.v);
      Float_t tmpZcomp = (zPos - trk->vtx.z);
      Float_t tmpAngle = GetAngle(tmpVcomp,tmpZcomp);
      Int_t weighting = (Int_t) hitScor;
      for(Int_t i=0;i<weighting;i++){
	allHitSummedAnglesV+=tmpAngle;
	allHitVangleCount++;
      }
    }
    if(IsUPln(pln)){
      Float_t tmpUcomp = (tPos - trk->vtx.u);
      Float_t tmpZcomp = (zPos - trk->vtx.z);
      Float_t tmpAngle = GetAngle(tmpUcomp,tmpZcomp);
      Int_t weighting = (Int_t) hitScor;
      for(Int_t i=0;i<weighting;i++){
	allHitSummedAnglesU+=tmpAngle;
	allHitUangleCount++;
      }
    }

    if(((stp->ph0.pe)+(stp->ph1.pe))<1.5) continue;
    if(((stp->z)-(trk->vtx.z))>2) continue;

    if(!IsTrkHit(zPos,tPos,trk,mb)){

      numVHShits++;
      remainingPH+=hitScor;
      if(IsVPln(pln)) HTVpln.FillHough(zPos,tPos);
      if(IsUPln(pln)) HTUpln.FillHough(zPos,tPos); 

      if(IsVPln(pln)){
	Float_t tmpVcomp = (tPos - trk->vtx.v);
	Float_t tmpZcomp = (zPos - trk->vtx.z);
	Float_t tmpAngle = GetAngle(tmpVcomp,tmpZcomp);
	Float_t protRangeAlongHitDir = (cos(tmpAngle)*predProtRange);
	if(protRangeAlongHitDir<(pln-trackVtxPlane) && protRangeAlongHitDir>0) phInProtRange+=hitScor;
	numAnglesV++;
	if(numAnglesV==1){
	  lowAngleV = highAngleV = tmpAngle;
	}
	if(numAnglesV>1){
	  if(tmpAngle<lowAngleV) lowAngleV = tmpAngle;
	  if(tmpAngle>highAngleV) highAngleV = tmpAngle;
	}
	Int_t weighting = (Int_t) hitScor;
	for(Int_t i=0;i<weighting;i++){
	  VHShitSummedAnglesV+=tmpAngle;
	  VHShitVangleCount++;
	}
	if(tmpAngle>ev0_5sigLowV && tmpAngle<ev0_5sigHighV){
	  ph0_5sigV+=hitScor;
	  nh0_5sigV++;
	}
	else if(tmpAngle>ev1_0sigLowV && tmpAngle<ev1_0sigHighV){
	  ph1_0sigV+=hitScor;
	  nh1_0sigV++;
	}
	else if(tmpAngle>ev1_5sigLowV && tmpAngle<ev1_5sigHighV){
	  ph1_5sigV+=hitScor;
	  nh1_5sigV++;
	}
	else if(tmpAngle>ev2_0sigLowV && tmpAngle<ev2_0sigHighV){
	  ph2_0sigV+=hitScor;
	  nh2_0sigV++;
	}
      }

      if(IsUPln(pln)){
	Float_t tmpUcomp = (tPos - trk->vtx.u);
	Float_t tmpZcomp = (zPos - trk->vtx.z);
	Float_t tmpAngle = GetAngle(tmpUcomp,tmpZcomp);
	Float_t protRangeAlongHitDir = (cos(tmpAngle)*predProtRange);
	if(protRangeAlongHitDir<(pln-trackVtxPlane) && protRangeAlongHitDir>0) phInProtRange+=hitScor;
	numAnglesU++;
	if(numAnglesU==1){
	  lowAngleU = highAngleU = tmpAngle;
	}
	if(numAnglesU>1){
	  if(tmpAngle<lowAngleU) lowAngleU = tmpAngle;
	  if(tmpAngle>highAngleU) highAngleU = tmpAngle;
	}
	Int_t weighting = (Int_t) hitScor;
	for(Int_t i=0;i<weighting;i++){
	  VHShitSummedAnglesU+=tmpAngle;
	  VHShitUangleCount++;
	}
	if(tmpAngle>ev0_5sigLowU && tmpAngle<ev0_5sigHighU){
	  ph0_5sigU+=hitScor;
	  nh0_5sigU++;
	}
	else if(tmpAngle>ev1_0sigLowU && tmpAngle<ev1_0sigHighU){
	  ph1_0sigU+=hitScor;
	  nh1_0sigU++;
	}
	else if(tmpAngle>ev1_5sigLowU && tmpAngle<ev1_5sigHighU){
	  ph1_5sigU+=hitScor;
	  nh1_5sigU++;
	}
	else if(tmpAngle>ev2_0sigLowU && tmpAngle<ev2_0sigHighU){
	  ph2_0sigU+=hitScor;
	  nh2_0sigU++;
	}
      }
      
    }

    for(Int_t ishw=0;ishw<evt->nshower;ishw++){
      int shwindex = evt->shw[ishw];    
      const NtpSRShower* shwr = mb->GetShower(shwindex);
      if(!shwr) continue;
      if(IsShowHit(zPos,tPos,shwr,mb) && IsTrkHit(zPos,tPos,trk,mb)){

	if((hitScor-oneMIPinScor)>0.0) {
	  numVHShits++;
	  remainingPH+=(hitScor-oneMIPinScor);
	  if(IsVPln(pln)) HTVpln.FillHough(zPos,tPos);
	  if(IsUPln(pln)) HTUpln.FillHough(zPos,tPos); 

	  if(IsVPln(pln)){
	    Float_t tmpVcomp = (tPos - trk->vtx.v);
	    Float_t tmpZcomp = (zPos - trk->vtx.z);
	    Float_t tmpAngle = GetAngle(tmpVcomp,tmpZcomp);
	    Float_t protRangeAlongHitDir = (cos(tmpAngle)*predProtRange);
	    if(protRangeAlongHitDir<(pln-trackVtxPlane) && protRangeAlongHitDir>0) phInProtRange+=hitScor;
	    numAnglesV++;
	    if(numAnglesV==1){
	      lowAngleV = highAngleV = tmpAngle;
	    }
	    if(numAnglesV>1){
	      if(tmpAngle<lowAngleV) lowAngleV = tmpAngle;
	      if(tmpAngle>highAngleV) highAngleV = tmpAngle;
	    }
	    Int_t weighting = (Int_t) (hitScor-oneMIPinScor);
	    for(Int_t i=0;i<weighting;i++){
	      VHShitSummedAnglesV+=tmpAngle;
	      VHShitVangleCount++;
	    }
	    if(tmpAngle>ev0_5sigLowV && tmpAngle<ev0_5sigHighV){
	      ph0_5sigV+=(hitScor-oneMIPinScor);
	      nh0_5sigV++;
	    }
	    else if(tmpAngle>ev1_0sigLowV && tmpAngle<ev1_0sigHighV){
	      ph1_0sigV+=(hitScor-oneMIPinScor);
	      nh1_0sigV++;
	    }
	    else if(tmpAngle>ev1_5sigLowV && tmpAngle<ev1_5sigHighV){
	      ph1_5sigV+=(hitScor-oneMIPinScor);
	      nh1_5sigV++;
	    }
	    else if(tmpAngle>ev2_0sigLowV && tmpAngle<ev2_0sigHighV){
	      ph2_0sigV+=(hitScor-oneMIPinScor);
	      nh2_0sigV++;
	    }
	  }
	  
	  if(IsUPln(pln)){
	    Float_t tmpUcomp = (tPos - trk->vtx.u);
	    Float_t tmpZcomp = (zPos - trk->vtx.z);
	    Float_t tmpAngle = GetAngle(tmpUcomp,tmpZcomp);
	    Float_t protRangeAlongHitDir = (cos(tmpAngle)*predProtRange);
	    if(protRangeAlongHitDir<(pln-trackVtxPlane) && protRangeAlongHitDir>0) phInProtRange+=hitScor;
	    numAnglesU++;
	    if(numAnglesU==1){
	      lowAngleU = highAngleU = tmpAngle;
	    }
	    if(numAnglesU>1){
	      if(tmpAngle<lowAngleU) lowAngleU = tmpAngle;
	      if(tmpAngle>highAngleU) highAngleU = tmpAngle;
	    }
	    Int_t weighting = (Int_t) (hitScor-oneMIPinScor);
	    for(Int_t i=0;i<weighting;i++){
	      VHShitSummedAnglesU+=tmpAngle;
	      VHShitUangleCount++;
	    }
	    if(tmpAngle>ev0_5sigLowU && tmpAngle<ev0_5sigHighU){
	      ph0_5sigU+=(hitScor-oneMIPinScor);
	      nh0_5sigU++;
	    }
	    else if(tmpAngle>ev1_0sigLowU && tmpAngle<ev1_0sigHighU){
	      ph1_0sigU+=(hitScor-oneMIPinScor);
	      nh1_0sigU++;
	    }
	    else if(tmpAngle>ev1_5sigLowU && tmpAngle<ev1_5sigHighU){
	      ph1_5sigU+=(hitScor-oneMIPinScor);
	      nh1_5sigU++;
	    }
	    else if(tmpAngle>ev2_0sigLowU && tmpAngle<ev2_0sigHighU){
	      ph2_0sigU+=(hitScor-oneMIPinScor);
	      nh2_0sigU++;
	    }
	  }
	  
	}

      }
    }

  }

  houghPeak = (HTUpln.GetPeakHeight())+(HTVpln.GetPeakHeight());
  rePH = remainingPH;
  vhsPH = remainingPH;
  totPH = totalEventPH;
  if(totalEventPH>0.0) rePHfrac = (remainingPH/totalEventPH);
  else rePHfrac = 0.0;
  nVHShits = numVHShits;

  if(allHitVangleCount>0){
    allHitAngV = (allHitSummedAnglesV / (Float_t) allHitVangleCount);
  }
  if(allHitUangleCount>0){
    allHitAngU = (allHitSummedAnglesU / (Float_t) allHitUangleCount);
  }
  if(VHShitVangleCount>0){
    VHSangV = (VHShitSummedAnglesV / (Float_t) VHShitVangleCount);
  }
  if(VHShitUangleCount>0){
    VHSangU = (VHShitSummedAnglesU / (Float_t) VHShitUangleCount);
  }

  houghRmsU = (HTUpln.houghSpace)->GetRMS(1);
  houghRmsV = (HTVpln.houghSpace)->GetRMS(1);

  phIn0_5sigCone = (ph0_5sigU+ph0_5sigV);
  phIn1_0sigCone = (ph1_0sigU+ph1_0sigV);
  phIn1_5sigCone = (ph1_5sigU+ph1_5sigV);
  phIn2_0sigCone = (ph2_0sigU+ph2_0sigV);

  nhIn0_5sigCone = (nh0_5sigU+nh0_5sigV);
  nhIn1_0sigCone = (nh1_0sigU+nh1_0sigV);
  nhIn1_5sigCone = (nh1_5sigU+nh1_5sigV);
  nhIn2_0sigCone = (nh2_0sigU+nh2_0sigV);

  //userVar = 0.0;

  if(numAnglesU>1 && numAnglesV>1){
    shwAngSpread = (((highAngleV-lowAngleV)+(highAngleU-lowAngleU))/2.0);
  }
  if(numAnglesU>1 && numAnglesV<=1) shwAngSpread = (highAngleU-lowAngleU);
  if(numAnglesV>1 && numAnglesU<=1) shwAngSpread = (highAngleV-lowAngleV);

  phInProtR = phInProtRange;

  return true;

}

//********************************************************
Float_t NewMadQEID::GetAngle(Float_t compUV, Float_t compZ)
{

  Float_t cosTheta = -10.0;

  if((sqrt(compUV*compUV + compZ*compZ))!=0) {
    cosTheta = (compZ/(sqrt(compUV*compUV + compZ*compZ)));
  }
  else return -10.0;

  Int_t dirCase = GetDirectionCase(compUV,compZ);

  if(dirCase==-1) return -10.0;
  else if(dirCase==1){
    return (acos(cosTheta));
  }
  else if(dirCase==2){
    return (-acos(cosTheta));
  }
  else if(dirCase==3){
    return (-PI+acos(cosTheta)); 
  }
  else if(dirCase==4){
    return (PI-acos(cosTheta));
  }
  else if(dirCase==5){
    return 0.0;
  }
  else if(dirCase==6){
    return PI;
  }
  else if(dirCase==7){
    return (PI/2.0);
  }
  else if(dirCase==8){
    return (-PI/2.0);
  }
  else return -10.0;

}

//********************************************************
Int_t NewMadQEID::GetDirectionCase(Float_t compUV, Float_t compZ)
{

  if(compUV==0){
    if(compZ>0) return 5;
    if(compZ<0) return 6;
  }
  else {
    if(compUV>0){
      if(compZ>0) return 1;
      if(compZ<0) return 4;
    }
    if(compUV<0){
      if(compZ>0) return 2;
      if(compZ<0) return 3;
    }
  }

  if(compZ==0){
    if(compUV>0) return 7;
    if(compUV<0) return 8;
  }

  return -1;

}

//********************************************************
Float_t NewMadQEID::CalcQEID(MadBase* mb, Int_t evtindex, Float_t reco_enu, 
			     Float_t reco_qe_enu, Float_t reco_emu, Float_t recoW2,
                             Int_t ntotss)
{

  if(!pdfsNormed){
    NormalisePDFs();
  }

  if(reco_enu>=maxE) return -10.0;

  if(!CalcQEVars(mb,evtindex,reco_enu,reco_emu)){
    return -15.0;
  }

  recoWsqd = recoW2;
  nSubShows = ntotss;
  qeKinEnuDif = (reco_enu-reco_qe_enu);

  // declare events automatically QE or non-QE for certain variable values
  if(recoWsqd>6.0) return -5.0;

  // calculate compound probabilities and ID parameter
  Int_t iBin = GetEnergyBin(reco_enu);

  Float_t QEnShow = 0.0;
  Float_t QEhPeak = 0.0;
  Float_t QErecoW2 = 0.0;
  Float_t QErePHfrac = 0.0;
  Float_t QEnVHShits = 0.0;

  Float_t QEnSubShow = 0.0;
  Float_t QEvhsPH = 0.0;
  Float_t QEqeKinEnuDif = 0.0;
  Float_t QEuserVar = 0.0;

  Float_t BGnShow = 0.0;
  Float_t BGhPeak = 0.0;
  Float_t BGrecoW2 = 0.0;
  Float_t BGrePHfrac = 0.0;
  Float_t BGnVHShits = 0.0;

  Float_t BGnSubShow = 0.0;
  Float_t BGvhsPH = 0.0;
  Float_t BGqeKinEnuDif = 0.0;
  Float_t BGuserVar = 0.0;

  // QE probs
  if(useNshows) QEnShow = nShowQE[iBin]->GetBinContent(nShowQE[iBin]->FindBin(nShows));
  if(useHoughPeak) QEhPeak = hPeakQE[iBin]->GetBinContent(hPeakQE[iBin]->FindBin(houghPeak));
  if(useRecoWsqd) QErecoW2 = recoW2QE[iBin]->GetBinContent(recoW2QE[iBin]->FindBin(recoWsqd));
  if(useRePHfrac) QErePHfrac = rePHfracQE[iBin]->GetBinContent(rePHfracQE[iBin]->FindBin(rePHfrac));
  if(useNVHShits) QEnVHShits = nVHShitsQE[iBin]->GetBinContent(nVHShitsQE[iBin]->FindBin(nVHShits));

  if(useNsubshows) QEnSubShow = nSSQE[iBin]->GetBinContent(nSSQE[iBin]->FindBin(nSubShows));
  if(useVhsPH) QEvhsPH = vhsPHQE[iBin]->GetBinContent(vhsPHQE[iBin]->FindBin(vhsPH));
  if(useQeKinEnuDif) QEqeKinEnuDif = qeKinEnuDifQE[iBin]->GetBinContent(qeKinEnuDifQE[iBin]->FindBin(qeKinEnuDif));
  if(useUserVariable) QEuserVar = userVarQE[iBin]->GetBinContent(userVarQE[iBin]->FindBin(userVar));

  // BG probs
  if(useNshows) BGnShow = nShowBG[iBin]->GetBinContent(nShowBG[iBin]->FindBin(nShows));
  if(useHoughPeak) BGhPeak = hPeakBG[iBin]->GetBinContent(hPeakBG[iBin]->FindBin(houghPeak));
  if(useRecoWsqd) BGrecoW2 = recoW2BG[iBin]->GetBinContent(recoW2BG[iBin]->FindBin(recoWsqd));
  if(useRePHfrac) BGrePHfrac = rePHfracBG[iBin]->GetBinContent(rePHfracBG[iBin]->FindBin(rePHfrac));
  if(useNVHShits) BGnVHShits = nVHShitsBG[iBin]->GetBinContent(nVHShitsBG[iBin]->FindBin(nVHShits));

  if(useNsubshows) BGnSubShow = nSSBG[iBin]->GetBinContent(nSSBG[iBin]->FindBin(nSubShows));
  if(useVhsPH) BGvhsPH = vhsPHBG[iBin]->GetBinContent(vhsPHBG[iBin]->FindBin(vhsPH));
  if(useQeKinEnuDif) BGqeKinEnuDif = qeKinEnuDifBG[iBin]->GetBinContent(qeKinEnuDifBG[iBin]->FindBin(qeKinEnuDif));
  if(useUserVariable) BGuserVar = userVarBG[iBin]->GetBinContent(userVarBG[iBin]->FindBin(userVar));

  Int_t numQElikeVarsZero = 0;
  Bool_t isQEPHfracZero = false;
  Bool_t isQEw2zero = false;
  Bool_t isQEnshwZero = false;
  Bool_t isQEhpkZero = false;
  Bool_t isQEnHitsZero = false;
  Bool_t isQEnSubShowZero = false;
  Bool_t isQEvhsPHZero = false;
  Bool_t isQEqeKinEnuDifZero = false;
  Bool_t isQEuserVarZero = false;
  if(useRePHfrac && QErePHfrac==0){
    numQElikeVarsZero++;
    isQEPHfracZero = true;
  }
  if(useRecoWsqd && QErecoW2==0){
    numQElikeVarsZero++;
    isQEw2zero = true;
  }
  if(useNshows && QEnShow==0){
    numQElikeVarsZero++;
    isQEnshwZero = true;
  }
  if(useHoughPeak && QEhPeak==0){
    numQElikeVarsZero++;
    isQEhpkZero = true;
  }
  if(useNVHShits && QEnVHShits==0){
    numQElikeVarsZero++;
    isQEnHitsZero = true;
  }
  if(useNsubshows && QEnSubShow==0){
    numQElikeVarsZero++;
    isQEnSubShowZero = true;
  }
  if(useVhsPH && QEvhsPH==0){
    numQElikeVarsZero++;
    isQEvhsPHZero = true;
  }
  if(useQeKinEnuDif && QEqeKinEnuDif==0){
    numQElikeVarsZero++;
    isQEqeKinEnuDifZero = true;
  }
  if(useUserVariable && QEuserVar==0){
    numQElikeVarsZero++;
    isQEuserVarZero = true;
  }

  Int_t numBGlikeVarsZero = 0;
  Bool_t isBGPHfracZero = false;
  Bool_t isBGw2zero = false;
  Bool_t isBGnshwZero = false;
  Bool_t isBGhpkZero = false;
  Bool_t isBGnHitsZero = false;
  Bool_t isBGnSubShowZero = false;
  Bool_t isBGvhsPHZero = false;
  Bool_t isBGqeKinEnuDifZero = false;
  Bool_t isBGuserVarZero = false;
  if(useRePHfrac && BGrePHfrac==0){
    numBGlikeVarsZero++;
    isBGPHfracZero = true;
  }
  if(useRecoWsqd && BGrecoW2==0){
    numBGlikeVarsZero++;
    isBGw2zero = true;
  }
  if(useNshows && BGnShow==0){
    numBGlikeVarsZero++;
    isBGnshwZero = true;
  }
  if(useHoughPeak && BGhPeak==0){
    numBGlikeVarsZero++;
    isBGhpkZero = true;
  }
  if(useNVHShits && BGnVHShits==0){
    numBGlikeVarsZero++;
    isBGnHitsZero = true;
  }
  if(useNsubshows && BGnSubShow==0){
    numBGlikeVarsZero++;
    isBGnSubShowZero = true;
  }
  if(useVhsPH && BGvhsPH==0){
    numBGlikeVarsZero++;
    isBGvhsPHZero = true;
  }
  if(useQeKinEnuDif && BGqeKinEnuDif==0){
    numBGlikeVarsZero++;
    isBGqeKinEnuDifZero = true;
  }
  if(useUserVariable && BGuserVar==0){
    numBGlikeVarsZero++;
    isBGuserVarZero = true;
  }

  Float_t likliQE = 0.0;
  Float_t likliNotQE = 0.0;

  if(numQElikeVarsZero>=2) likliQE = 0.0;
  if(numBGlikeVarsZero>=2) likliNotQE = 0.0;

  if(numQElikeVarsZero==1){
    if(isQEPHfracZero) {
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEw2zero) {
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEnshwZero) {
      likliQE = 1.0;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEhpkZero) {
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEnHitsZero) {
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEnSubShowZero){
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEvhsPHZero){
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEqeKinEnuDifZero){
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useUserVariable) likliQE*=QEuserVar;
    }
    if(isQEuserVarZero){
      likliQE = 1.0;
      if(useNshows) likliQE*=QEnShow;
      if(useHoughPeak) likliQE*=QEhPeak;
      if(useRecoWsqd) likliQE*=QErecoW2;
      if(useRePHfrac) likliQE*=QErePHfrac;
      if(useNVHShits) likliQE*=QEnVHShits;
      if(useNsubshows) likliQE*=QEnSubShow;
      if(useVhsPH) likliQE*=QEvhsPH;
      if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
    }

  }
  if(numBGlikeVarsZero==1){
    if(isBGPHfracZero) {
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGw2zero) {
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGnshwZero) {
      likliNotQE = 1.0;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGhpkZero) {
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGnHitsZero) {
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGnSubShowZero){
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGvhsPHZero){
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGqeKinEnuDifZero){
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useUserVariable) likliNotQE*=BGuserVar;
    }
    if(isBGuserVarZero){
      likliNotQE = 1.0;
      if(useNshows) likliNotQE*=BGnShow;
      if(useHoughPeak) likliNotQE*=BGhPeak;
      if(useRecoWsqd) likliNotQE*=BGrecoW2;
      if(useRePHfrac) likliNotQE*=BGrePHfrac;
      if(useNVHShits) likliNotQE*=BGnVHShits;
      if(useNsubshows) likliNotQE*=BGnSubShow;
      if(useVhsPH) likliNotQE*=BGvhsPH;
      if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
    }

  }

  if(numQElikeVarsZero==0){
    likliQE = 1.0;
    if(useNshows) likliQE*=QEnShow;
    if(useHoughPeak) likliQE*=QEhPeak;
    if(useRecoWsqd) likliQE*=QErecoW2;
    if(useRePHfrac) likliQE*=QErePHfrac;
    if(useNVHShits) likliQE*=QEnVHShits;
    if(useNsubshows) likliQE*=QEnSubShow;
    if(useVhsPH) likliQE*=QEvhsPH;
    if(useQeKinEnuDif) likliQE*=QEqeKinEnuDif;
    if(useUserVariable) likliQE*=QEuserVar;
  }
  if(numBGlikeVarsZero==0){
    likliNotQE = 1.0;
    if(useNshows) likliNotQE*=BGnShow;
    if(useHoughPeak) likliNotQE*=BGhPeak;
    if(useRecoWsqd) likliNotQE*=BGrecoW2;
    if(useRePHfrac) likliNotQE*=BGrePHfrac;
    if(useNVHShits) likliNotQE*=BGnVHShits;
    if(useNsubshows) likliNotQE*=BGnSubShow;
    if(useVhsPH) likliNotQE*=BGvhsPH;
    if(useQeKinEnuDif) likliNotQE*=BGqeKinEnuDif;
    if(useUserVariable) likliNotQE*=BGuserVar;
  }

  if(likliQE==0 && likliNotQE==0) return -10.0;
  if(likliQE==0 && likliNotQE>0) return -5.0;
  if(likliQE>0 && likliNotQE==0) return 5.0;
  Float_t EID = -(sqrt(-log(likliQE)))+(sqrt(-log(likliNotQE)));
  return EID;  

}

//********************************************************
void NewMadQEID::NormalisePDFs()
{

  if(pdfsNormed) return;

  Float_t integ;

  for(Int_t i=0;i<30;i++){

    if(useRePHfrac){
      integ=rePHfracQE[i]->Integral(0,((rePHfracQE[i]->GetNbinsX())+1));
      if(integ!=0) rePHfracQE[i]->Scale(1./integ);
      integ=rePHfracBG[i]->Integral(0,((rePHfracBG[i]->GetNbinsX())+1));
      if(integ!=0) rePHfracBG[i]->Scale(1./integ);
    }

    if(useRecoWsqd){
      integ=recoW2QE[i]->Integral(0,((recoW2QE[i]->GetNbinsX())+1));
      if(integ!=0) recoW2QE[i]->Scale(1./integ);
      integ=recoW2BG[i]->Integral(0,((recoW2BG[i]->GetNbinsX())+1));
      if(integ!=0) recoW2BG[i]->Scale(1./integ);
    }

    if(useNshows){
      integ=nShowQE[i]->Integral(0,((nShowQE[i]->GetNbinsX())+1));
      if(integ!=0) nShowQE[i]->Scale(1./integ);
      integ=nShowBG[i]->Integral(0,((nShowBG[i]->GetNbinsX())+1));
      if(integ!=0) nShowBG[i]->Scale(1./integ);
    }

    if(useHoughPeak){
      integ=hPeakQE[i]->Integral(0,((hPeakQE[i]->GetNbinsX())+1));
      if(integ!=0) hPeakQE[i]->Scale(1./integ);
      integ=hPeakBG[i]->Integral(0,((hPeakBG[i]->GetNbinsX())+1));
      if(integ!=0) hPeakBG[i]->Scale(1./integ);
    }

    if(useNVHShits){
      integ=nVHShitsQE[i]->Integral(0,((nVHShitsQE[i]->GetNbinsX())+1));
      if(integ!=0) nVHShitsQE[i]->Scale(1./integ);
      integ=nVHShitsBG[i]->Integral(0,((nVHShitsBG[i]->GetNbinsX())+1));
      if(integ!=0) nVHShitsBG[i]->Scale(1./integ);
    }

    if(useNsubshows){
      integ=nSSQE[i]->Integral(0,((nSSQE[i]->GetNbinsX())+1));
      if(integ!=0) nSSQE[i]->Scale(1./integ);
      integ=nSSBG[i]->Integral(0,((nSSBG[i]->GetNbinsX())+1));
      if(integ!=0) nSSBG[i]->Scale(1./integ);
    }

    if(useVhsPH){
      integ=vhsPHQE[i]->Integral(0,((vhsPHQE[i]->GetNbinsX())+1));
      if(integ!=0) vhsPHQE[i]->Scale(1./integ);
      integ=vhsPHBG[i]->Integral(0,((vhsPHBG[i]->GetNbinsX())+1));
      if(integ!=0) vhsPHBG[i]->Scale(1./integ);
    }

    if(useQeKinEnuDif){
      integ=qeKinEnuDifQE[i]->Integral(0,((qeKinEnuDifQE[i]->GetNbinsX())+1));
      if(integ!=0) qeKinEnuDifQE[i]->Scale(1./integ);
      integ=qeKinEnuDifBG[i]->Integral(0,((qeKinEnuDifBG[i]->GetNbinsX())+1));
      if(integ!=0) qeKinEnuDifBG[i]->Scale(1./integ);
    }

    if(useUserVariable){
      integ=userVarQE[i]->Integral(0,((userVarQE[i]->GetNbinsX())+1));
      if(integ!=0) userVarQE[i]->Scale(1./integ);
      integ=userVarBG[i]->Integral(0,((userVarBG[i]->GetNbinsX())+1));
      if(integ!=0) userVarBG[i]->Scale(1./integ);
    }

  }

  pdfsNormed=true;

  return;

}

//********************************************************
void NewMadQEID::FillPDFs(Int_t cc_nc, Int_t process, Int_t inu, MadBase* mb, 
			  Int_t evtindex, Float_t reco_enu, Float_t reco_qe_enu, 
			  Float_t reco_emu, Float_t recoW2, Int_t ntotss)
{

  static Bool_t firstFill = true;
  if(firstFill){
    InitPDFs();
    firstFill = false;
  }

  if(reco_enu>=maxE) return;
  if(!CalcQEVars(mb,evtindex,reco_enu,reco_emu)) return;
  recoWsqd = recoW2;
  nSubShows = ntotss;
  qeKinEnuDif = (reco_enu-reco_qe_enu);

  Int_t iBin = GetEnergyBin(reco_enu);

  if(cc_nc==0 || inu==-14){

    if(useNshows) nShowBG[iBin]->Fill(nShows);
    if(useHoughPeak) hPeakBG[iBin]->Fill(houghPeak);
    if(useRecoWsqd) recoW2BG[iBin]->Fill(recoWsqd);
    if(useRePHfrac) rePHfracBG[iBin]->Fill(rePHfrac);
    if(useNVHShits) nVHShitsBG[iBin]->Fill(nVHShits);
    if(useNsubshows) nSSBG[iBin]->Fill(nSubShows);
    if(useVhsPH) vhsPHBG[iBin]->Fill(vhsPH);
    if(useQeKinEnuDif) qeKinEnuDifBG[iBin]->Fill(qeKinEnuDif);
    if(useUserVariable) userVarBG[iBin]->Fill(userVar);

  }

  if(cc_nc==1 && process!=1001 && inu==14){

    if(useNshows) nShowBG[iBin]->Fill(nShows);
    if(useHoughPeak) hPeakBG[iBin]->Fill(houghPeak);
    if(useRecoWsqd) recoW2BG[iBin]->Fill(recoWsqd);
    if(useRePHfrac) rePHfracBG[iBin]->Fill(rePHfrac);
    if(useNVHShits) nVHShitsBG[iBin]->Fill(nVHShits);
    if(useNsubshows) nSSBG[iBin]->Fill(nSubShows);
    if(useVhsPH) vhsPHBG[iBin]->Fill(vhsPH);
    if(useQeKinEnuDif) qeKinEnuDifBG[iBin]->Fill(qeKinEnuDif);
    if(useUserVariable) userVarBG[iBin]->Fill(userVar);

  }

  if(cc_nc==1 && process==1001 && inu==14){

    if(useNshows) nShowQE[iBin]->Fill(nShows);
    if(useHoughPeak) hPeakQE[iBin]->Fill(houghPeak);
    if(useRecoWsqd) recoW2QE[iBin]->Fill(recoWsqd);
    if(useRePHfrac) rePHfracQE[iBin]->Fill(rePHfrac);
    if(useNVHShits) nVHShitsQE[iBin]->Fill(nVHShits);
    if(useNsubshows) nSSQE[iBin]->Fill(nSubShows);
    if(useVhsPH) vhsPHQE[iBin]->Fill(vhsPH);
    if(useQeKinEnuDif) qeKinEnuDifQE[iBin]->Fill(qeKinEnuDif);
    if(useUserVariable) userVarQE[iBin]->Fill(userVar);

  }

  return;

}

//********************************************************
void NewMadQEID::WritePDFs(const char* name)
{

  TDirectory* dir = gDirectory;
  TFile* f = new TFile(name,"recreate");
  dir->cd();

  for(Int_t i=0;i<30;i++){

    if(useRePHfrac){
      rePHfracQE[i]->SetDirectory(f);
      rePHfracBG[i]->SetDirectory(f);
    }

    if(useRecoWsqd){
      recoW2QE[i]->SetDirectory(f);
      recoW2BG[i]->SetDirectory(f);
    }

    if(useNshows){
      nShowQE[i]->SetDirectory(f);
      nShowBG[i]->SetDirectory(f);
    }

    if(useHoughPeak){
      hPeakQE[i]->SetDirectory(f);
      hPeakBG[i]->SetDirectory(f);
    }

    if(useNVHShits){
      nVHShitsQE[i]->SetDirectory(f);
      nVHShitsBG[i]->SetDirectory(f);
    }

    if(useNsubshows){
      nSSQE[i]->SetDirectory(f);
      nSSBG[i]->SetDirectory(f);
    }

    if(useVhsPH){
      vhsPHQE[i]->SetDirectory(f);
      vhsPHBG[i]->SetDirectory(f);
    }

    if(useQeKinEnuDif){
      qeKinEnuDifQE[i]->SetDirectory(f);
      qeKinEnuDifBG[i]->SetDirectory(f);
    }

    if(useUserVariable){
      userVarQE[i]->SetDirectory(f);
      userVarBG[i]->SetDirectory(f);
    }

    if(trainingHists[i]) trainingHists[i]->SetDirectory(f);

  } 

  if(cutsHist) cutsHist->SetDirectory(f);
  nBinsHist->SetDirectory(f);
  eBins->SetDirectory(f);
  f->Write();

}

//********************************************************
void NewMadQEID::ReadPDFs(const char* name)
{

  TFile* QEfile = new TFile(name,"READ");

  if(QEfile->IsOpen() && !QEfile->IsZombie()){

    std::cout<<"Read MadQEID PDFs from : "<<name<<std::endl;
    QEfile->cd();

    for(Int_t i=0;i<30;i++){
      char nam[30];

      sprintf(nam,"rePHfracQE[%d]",i);
      if(useRePHfrac) rePHfracQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"rePHfracBG[%d]",i);
      if(useRePHfrac) rePHfracBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"recoW2QE[%d]",i);
      if(useRecoWsqd) recoW2QE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"recoW2BG[%d]",i);
      if(useRecoWsqd) recoW2BG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"hPeakQE[%d]",i);
      if(useHoughPeak) hPeakQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"hPeakBG[%d]",i);
      if(useHoughPeak) hPeakBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"nShowQE[%d]",i);
      if(useNshows) nShowQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"nShowBG[%d]",i);
      if(useNshows) nShowBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"nVHShitsQE[%d]",i);
      if(useNVHShits) nVHShitsQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"nVHShitsBG[%d]",i);
      if(useNVHShits) nVHShitsBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"nSSQE[%d]",i);
      if(useNsubshows) nSSQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"nSSBG[%d]",i);
      if(useNsubshows) nSSBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"vhsPHQE[%d]",i);
      if(useVhsPH) vhsPHQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"vhsPHBG[%d]",i);
      if(useVhsPH) vhsPHBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"qeKinEnuDifQE[%d]",i);
      if(useQeKinEnuDif) qeKinEnuDifQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"qeKinEnuDifBG[%d]",i);
      if(useQeKinEnuDif) qeKinEnuDifBG[i] = (TH1F*) QEfile->Get(nam);
      
      sprintf(nam,"userVarQE[%d]",i);
      if(useUserVariable) userVarQE[i] = (TH1F*) QEfile->Get(nam);
      sprintf(nam,"userVarBG[%d]",i);
      if(useUserVariable) userVarBG[i] = (TH1F*) QEfile->Get(nam);
      
    }
  }
  else {
    std::cout<<"Failed to read MadQEID PDFs from : "<<name<<std::endl;
  }

  if(!pdfsNormed) NormalisePDFs();

  if((cutsHist = (TH1F*) QEfile->Get("cuthist"))) {};

  //cout << "-------------------------------------" << endl;
  //cout << "Calculating array of axis values..." << endl;
  TH1F* tmpeBins = (TH1F*) QEfile->Get("eBins");
  Float_t axArray[31];
  Int_t axArrayIndex = 0;
  for(Int_t i=0;i<tmpeBins->GetNbinsX();i++){
    if((tmpeBins->GetBinContent(i))!=0){
      if(axArrayIndex==31) continue;
      axArray[axArrayIndex] = tmpeBins->GetBinLowEdge(i);
      //cout << "Bin " << axArrayIndex  << " has low edge: " << axArray[axArrayIndex] << endl;
      axArrayIndex++;
    }
  }
  //cout << "End of true bins (+1) low edges." << endl;
  Int_t finalTrueBinIndex = axArrayIndex;
  for(Int_t i=finalTrueBinIndex;i<31;i++){
    axArray[i] = 0.0;
    //cout << "Non-Bin " << i << " has low edge: " << axArray[i] << endl;
  }
  delete tmpeBins;
  //cout << "------------------------------------" << endl;
  SetEnergyBinning(axArray);

}

//********************************************************
void NewMadQEID::FillTrainingHists(MadBase* mb, Int_t evtindex, Float_t reco_enu, 
				   Float_t reco_qe_enu, Float_t reco_emu, Float_t recoW2,
				   Int_t ntotss)
{

  static Bool_t firstTrainingFill = true;
  if(firstTrainingFill) {
    InitTrainingHists();
    firstTrainingFill = false;
  }

  if(reco_enu>=maxE) return;

  Int_t iBin = GetEnergyBin(reco_enu);
  Float_t iID = CalcQEID(mb,evtindex,reco_enu,reco_qe_enu,reco_emu,recoW2,ntotss);

  trainingHists[iBin]->Fill(iID);

  return;

}

//********************************************************
void NewMadQEID::MakeCutHist(Float_t eff)
{

  for(Int_t i=0;i<nBins;i++){

    Int_t totNumQE = (Int_t) trainingHists[i]->GetEntries();
    Float_t iBinStart = GetBinStart(i);
    Float_t iBinRange = GetBinRange(i);
    Float_t lookUpVal = (iBinStart + (iBinRange/2.0));

    // fill default cut value
    if(totNumQE==0){
      cutsHist->Fill(lookUpVal,0.5);
      continue;
    }

    else {

      Int_t currentNumQE = 0;
      Float_t currentEff = 0.0;
      Int_t nHistBins = trainingHists[i]->GetNbinsX();
      Int_t currentHistBin = nHistBins;

      while(currentEff<eff) {

	currentNumQE+=(Int_t) trainingHists[i]->GetBinContent(currentHistBin);
	if(currentNumQE==0) currentEff = 0.0;
	else {
	  currentEff = ((Float_t) currentNumQE / (Float_t) totNumQE);
	}
	currentHistBin--;
	if(currentHistBin==0) break;

      }

      if(currentHistBin==0) cutsHist->Fill(lookUpVal,4.5);  // still includes definite QE events
      else cutsHist->Fill(lookUpVal,trainingHists[i]->GetBinCenter(currentHistBin));

    }

  }

  return;

}

//********************************************************
Int_t NewMadQEID::PassMedQECut(Float_t reco_enu, Float_t QEIDval)
{
  
  if(reco_enu>=maxE) return -1;

  if(QEIDval==-10.0) return -1;
  else if(QEIDval==-5.0) return 0;
  else if(QEIDval==5.0) return 1;
  else {
    Int_t iBin = GetEnergyBin(reco_enu);
    Float_t iBinStart = GetBinStart(iBin);
    Float_t iBinRange = GetBinRange(iBin);
    Float_t lookUpVal = (iBinStart + (iBinRange/2.0));
    Float_t cutVal = cutsHist->GetBinContent(cutsHist->FindBin(lookUpVal));
    if(QEIDval>cutVal) return 1;
    else if(QEIDval<=cutVal) return 0;
  }

  return -1;

}

//********************************************************
Int_t NewMadQEID::GetEnergyBin(Float_t enu)
{

  if(enu<0.0) return -1;
  else if(enu>=maxE) return -1;
  else {
    Int_t tmpNbins = nBins;
    while(GetBinStart(tmpNbins)>enu){
      tmpNbins--;
    }
    return tmpNbins;
  }

}

//********************************************************
Float_t NewMadQEID::GetBinStart(Int_t bin)
{

  if(bin<0) return -1.0;
  else if(bin>nBins) return -1.0;
  else return binArray[bin];  // 'bin' index starts at zero

}

//********************************************************
Float_t NewMadQEID::GetBinEnd(Int_t bin)
{

  if(bin<0) return -1.0;
  else if(bin>nBins) return -1.0;
  else return binArray[bin+1];

}

//********************************************************
Float_t NewMadQEID::GetBinRange(Int_t bin)
{

  if(bin<0) return -1.0;
  else if(bin>nBins) return -1.0;
  else return (binArray[bin+1]-binArray[bin]);

}

//********************************************************
Float_t NewMadQEID::GetRecoWsqd()
{
  return recoWsqd;
}

//********************************************************
Int_t NewMadQEID::GetNshows()
{
  return nShows;
}

//********************************************************
Int_t NewMadQEID::GetHoughPeak()
{
  return houghPeak;
}

//********************************************************
Float_t NewMadQEID::GetRePHfrac()
{
  return rePHfrac;
}

//********************************************************
Int_t NewMadQEID::GetNVHShits()
{
  return nVHShits;
}

//********************************************************
Int_t NewMadQEID::GetNsubshows()
{
  return nSubShows;
}

//********************************************************
Float_t NewMadQEID::GetVHSPH(){
  return vhsPH;
}

//********************************************************
Float_t NewMadQEID::GetQEkinEnuDif(){
  return qeKinEnuDif;
}

//********************************************************
Float_t NewMadQEID::GetUserVar(){
  return userVar;
}

//********************************************************
Float_t  NewMadQEID::GetProtonAngU() 
{
  return protonAngU;
}

//********************************************************
Float_t  NewMadQEID::GetProtonAngV() 
{
  return protonAngV;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtAngU() 
{
  return stdHepProtAngU;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtAngV() 
{
  return stdHepProtAngV;
}

//********************************************************
Float_t  NewMadQEID::GetAllHitAngU() 
{
  return allHitAngU;
}

//********************************************************
Float_t  NewMadQEID::GetAllHitAngV() 
{
  return allHitAngV;
}

//********************************************************
Float_t  NewMadQEID::GetVHSangU() 
{
  return VHSangU;
}

//********************************************************
Float_t  NewMadQEID::GetVHSangV() 
{
  return VHSangV;
}

//********************************************************
Float_t  NewMadQEID::GetProtonUcomp() 
{
  return protonUcomp;
}

//********************************************************
Float_t  NewMadQEID::GetProtonVcomp() 
{
  return protonVcomp;
}

//********************************************************
Float_t  NewMadQEID::GetProtonZcomp() 
{
  return protonZcomp;
}

//********************************************************
Float_t  NewMadQEID::GetProtonEfromEvent() 
{
  return protonEfromEvent;
}

//********************************************************
Float_t  NewMadQEID::GetProtonEfromMassMom() 
{
  return protonEfromMassMom;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtUcomp() 
{
  return stdHepProtUcomp;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtVcomp() 
{
  return stdHepProtVcomp;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtZcomp() 
{
  return stdHepProtZcomp;
}

//********************************************************
Float_t  NewMadQEID::GetStdHepProtE() 
{
  return stdHepProtE;
}

//********************************************************
Float_t NewMadQEID::GetMuonUcomp()
{
  return muonUcomp;
}

//********************************************************
Float_t NewMadQEID::GetMuonVcomp()
{
  return muonVcomp;
}

//********************************************************
Float_t NewMadQEID::GetMuonZcomp()
{
  return muonZcomp;
}

//********************************************************
Float_t NewMadQEID::GetMuonE()
{
  return muonE;
}

//********************************************************
Float_t NewMadQEID::GetStdHepMuUcomp()
{
  return stdHepMuUcomp;
}

//********************************************************
Float_t NewMadQEID::GetStdHepMuVcomp()
{
  return stdHepMuVcomp;
}

//********************************************************
Float_t NewMadQEID::GetStdHepMuZcomp()
{
  return stdHepMuZcomp;
}

//********************************************************
Float_t NewMadQEID::GetStdHepMuE()
{
  return stdHepMuE;
}

//********************************************************
Float_t NewMadQEID::GetHoughRmsU()
{
  return houghRmsU;
}

//********************************************************
Float_t NewMadQEID::GetHoughRmsV()
{
  return houghRmsV;
}

//********************************************************
Float_t NewMadQEID::GetPhIn0_5sigCone()
{
  return phIn0_5sigCone;
}

//********************************************************
Float_t NewMadQEID::GetPhIn1_0sigCone()
{
  return phIn1_0sigCone;
}

//********************************************************
Float_t NewMadQEID::GetPhIn1_5sigCone()
{
  return phIn1_5sigCone;
}

//********************************************************
Float_t NewMadQEID::GetPhIn2_0sigCone()
{
  return phIn2_0sigCone;
}

//********************************************************
Float_t NewMadQEID::GetTeIn0_5sigCone()
{
  return teIn0_5sigCone;
}

//********************************************************
Float_t NewMadQEID::GetTeIn1_0sigCone()
{
  return teIn1_0sigCone;
}

//********************************************************
Float_t NewMadQEID::GetTeIn1_5sigCone()
{
  return teIn1_5sigCone;
}

//********************************************************
Float_t NewMadQEID::GetTeIn2_0sigCone()
{
  return teIn2_0sigCone;
}

//********************************************************
Float_t NewMadQEID::GetTeTot()
{
  return teTot;
}

//********************************************************
Int_t NewMadQEID::GetNhIn0_5sigCone()
{
  return nhIn0_5sigCone;
}

//********************************************************
Int_t NewMadQEID::GetNhIn1_0sigCone()
{
  return nhIn1_0sigCone;
}

//********************************************************
Int_t NewMadQEID::GetNhIn1_5sigCone()
{
  return nhIn1_5sigCone;
}

//********************************************************
Int_t NewMadQEID::GetNhIn2_0sigCone()
{
  return nhIn2_0sigCone;
}

//********************************************************
Float_t NewMadQEID::GetEvRePH()
{
  return rePH;
}

//********************************************************
Float_t NewMadQEID::GetEvTotPH()
{
  return totPH;
}

//********************************************************
Float_t NewMadQEID::GetUserVariable()
{
  return userVar;
}

//********************************************************
Float_t NewMadQEID::GetPhInProtRange()
{
  return phInProtR;
}

//********************************************************
Float_t NewMadQEID::GetShwAngSpread()
{
  return shwAngSpread;
}

//********************************************************
void NewMadQEID::UseNshows(Bool_t tof)
{
  useNshows = tof;
  return;
}

//********************************************************
void NewMadQEID::UseHoughPeak(Bool_t tof)
{
  useHoughPeak = tof;
  return;
}

//********************************************************
void NewMadQEID::UseRecoWsqd(Bool_t tof)
{
  useRecoWsqd = tof;
  return;
}

//********************************************************
void NewMadQEID::UseRePHfrac(Bool_t tof)
{
  useRePHfrac = tof;
  return;
}

//********************************************************
void NewMadQEID::UseNVHShits(Bool_t tof)
{
  useNVHShits = tof;
  return;
}

//********************************************************
void NewMadQEID::UseNsubshows(Bool_t tof)
{
  useNsubshows = tof;
  return;
}

//********************************************************
void NewMadQEID::UseVHSPH(Bool_t tof)
{
  useVhsPH = tof;
  return;
}

//********************************************************
void NewMadQEID::UseQEkinEnuDif(Bool_t tof)
{
  useQeKinEnuDif = tof;
  return;
}

//********************************************************
void NewMadQEID::UseUserVar(Bool_t tof)
{
  useUserVariable = tof;
  return;
}

//********************************************************
Bool_t NewMadQEID::IsTrkHit(Float_t zPos,Float_t tPos, const NtpSRTrack* ntpTrk, 
			    MadBase* madb)
{

  for(Int_t i=0;i<ntpTrk->nstrip;i++){
    const NtpSRStrip* strp = madb->GetStrip(ntpTrk->stp[i]);
    if(!strp) continue;
    if(strp->tpos==tPos && strp->z==zPos){
      return true;
    }
  }
  return false;

}

//********************************************************
Bool_t NewMadQEID::IsShowHit(Float_t zPos,Float_t tPos, const NtpSRShower* ntpShow, 
			     MadBase* madb)
{

  for(Int_t i=0;i<ntpShow->nstrip;i++){
    const NtpSRStrip* strp = madb->GetStrip(ntpShow->stp[i]);
    if(!strp) continue;
    if(strp->tpos==tPos && strp->z==zPos){
      return true;
    }
  }
  return false;

}

//********************************************************
Bool_t NewMadQEID::IsVPln(Int_t plane)
{

  if((plane%2)==0) return true;
  return false;

}

//********************************************************
Bool_t NewMadQEID::IsUPln(Int_t plane)
{

  if((plane%2)!=0) return true;
  return false;

}

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

#endif