VHS.cxx

Go to the documentation of this file.

// $Id: VHS.cxx,v 1.2 2007/04/27 16:40:45 arms Exp $
//
// VHS : Namespace containing functions useful for training
//       and analysis of event images
//
// arms@physics.umn.edu
#include "AnalysisNtuples/Module/VHS.h"
#include "AnalysisNtuples/Module/VHSevent.h"

#include "MinosObjectMap/MomNavigator.h"
#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "MCNtuple/NtpMCRecord.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "Conventions/PlaneView.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCTruth.h"
#include "TruthHelperNtuple/NtpTHEvent.h"

#include "TBenchmark.h"
#include "TButton.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH2D.h"
#include "TList.h"
#include "TMath.h"
#include "TPaveStats.h"
#include "TPaveText.h"
#include "TSystem.h"
#include "TTree.h"

#include "Riostream.h"
#include <fstream>
#include <iomanip>
#include <sstream>

//......................................................................

bool VHS::DrawEvent( int nEntry,  const char* treeName,
                     int nPlanes, int         nStrips )
{
  // Draw an event display for the given event from a TTree of VHSEvent

  TTree* theTree = (TTree*)gROOT->FindObject(treeName);

  if (!theTree) return false;
  const int kNumEntries = theTree->GetEntriesFast();

  if ( nEntry > kNumEntries ) return false;

  TCanvas* evtDisp = 0;

  if ( !gROOT->FindObject("VHSEvtDisp") ) {
    evtDisp = new TCanvas("VHSEvtDisp","VHS Event Display",0,0,800,600);
    evtDisp->Divide(2,2);
    evtDisp->GetPad(1)->SetPad(0.00,0.5,0.33,1.0);
    evtDisp->GetPad(2)->SetPad(0.33,0.5,1.00,1.0);
    evtDisp->GetPad(3)->SetPad(0.00,0.0,0.33,0.5);
    evtDisp->GetPad(4)->SetPad(0.33,0.0,1.00,0.5);

  } else {
    evtDisp = (TCanvas*)gROOT->FindObject("VHSEvtDisp");
    evtDisp->cd(1);
  }
  evtDisp->SetEditable(kTRUE);

  std::string nextCom = (nEntry+1 < kNumEntries)?
    std::string("VHS::DrawEvent(") + VHS::ToString(nEntry+1) + std::string(",\"") +
    std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
    VHS::ToString(nStrips) + std::string(");")                                           :
    std::string("VHS::DrawEvent(") + VHS::ToString(kNumEntries-1) + std::string(",\"") +
    std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
    VHS::ToString(nStrips) + std::string(");");

  std::string prevCom = (nEntry > 0)?
    std::string("VHS::DrawEvent(") + VHS::ToString(nEntry-1) + std::string(",\"") +
    std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
    VHS::ToString(nStrips) + std::string(");")                                           :
    std::string("VHS::DrawEvent(0,\"") +
    std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
    VHS::ToString(nStrips) + std::string(");");

  // Draw navigation buttons
  evtDisp->cd(3);
  evtDisp->GetPad(3)->GetListOfPrimitives()->Delete();
  evtDisp->GetPad(3)->Range(0,0,1,1);

  // Next button
  TButton* butNext =
    new TButton("Next Event",nextCom.c_str(),0.20,0.60,0.80,0.80);
  butNext->SetTextSize(0.3);

  // Previous button
  TButton* butPrev =
    new TButton("Prev Event",prevCom.c_str(),0.20,0.40,0.80,0.60);
  butPrev->SetTextSize(0.3);
    
  // Quit button
  TButton* butQuit = new TButton("Quit",".q",0.20,0.20,0.80,0.40);
  butQuit->SetTextSize(0.3);

  butNext->Draw();
  butPrev->Draw();
  butQuit->Draw();

  // Grab a copy of the event
  VHSevent* evt = new VHSevent();
  theTree->SetBranchAddress("VHSevent",&evt);
  theTree->GetEntry(nEntry);

  // Draw UZ view
  TH2D* uzImage = (TH2D*)gROOT->FindObject("UZImage");
  TH2D* vzImage = (TH2D*)gROOT->FindObject("VZImage");

  if ( uzImage ) { delete uzImage; uzImage = 0; }
  if ( vzImage ) { delete vzImage; vzImage = 0; }

  VHS::ToTH2D( VHS::FullVec( evt->stp, evt->ind, nPlanes, nStrips ),
               nPlanes, nStrips, "ZImage", "Z",
               uzImage, vzImage);

  evtDisp->cd(2);
  uzImage->SetXTitle("Plane");
  uzImage->SetYTitle("Strip");
  uzImage->SetStats(kFALSE);
  uzImage->Draw("COLZ");

  // Draw VZ view
  evtDisp->cd(4);
  vzImage->SetXTitle("Plane");
  vzImage->SetYTitle("Strip");
  vzImage->SetStats(kFALSE);
  vzImage->Draw("COLZ");

  // Draw Event information
  evtDisp->cd(1);
  evtDisp->GetPad(1)->GetListOfPrimitives()->Delete();

  bool   MC        = evt->MC;
  double MCenu     = evt->MCenu;
  int    MCiaction = evt->MCiaction;
  int    MCinu     = TMath::Abs(evt->MCinu);
  //int  MCrescode = evt->MCrescode;

  double dNC       = evt->dNC;
  double dCCe      = evt->dCCe;
  double dCCmu     = evt->dCCmu;
  double dCCtau    = evt->dCCtau;

  double pNC       = evt->pNC;
  double pCCe      = evt->pCCe;
  double pCCmu     = evt->pCCmu;
  double pCCtau    = evt->pCCtau;

  double GeV       = evt->GeV;
  double Sigcor    = evt->Sigcor;

  int    NumPln    = evt->NumPln;

  VHS::evtType evtt = VHS::GetEvtType( MCinu, MCiaction );

  // --> text boxes to dump info
  TPaveText* ptL = new TPaveText(0.02,0.02,0.5,0.45,"brNDC");
  ptL->SetLineWidth(2);
  ptL->SetBorderSize(1);
  ptL->SetTextSize(0.05);
  ptL->SetTextAlign(12);

  std::string li1L = std::string("dCCe = ")    + VHS::ToString(dCCe,7,5);
  std::string li2L = std::string("dCC#mu = ")  + VHS::ToString(dCCmu,7,5);
  std::string li3L = std::string("dCC#tau = ") + VHS::ToString(dCCtau,7,5);
  std::string li4L = std::string("dNC  = ")    + VHS::ToString(dNC,7,5);
  ptL->AddText(li1L.c_str());
  ptL->AddText(li2L.c_str());
  ptL->AddText(li3L.c_str());
  ptL->AddText(li4L.c_str());
  ptL->Draw();

  TPaveText* ptR = new TPaveText(0.5,0.02,0.98,0.45,"brNDC");
  ptR->SetLineWidth(2);
  ptR->SetBorderSize(1);
  ptR->SetTextSize(0.05);
  ptR->SetTextAlign(12);

  std::string li1R = std::string("pCCe = ")    + VHS::ToString(pCCe,7,5);
  std::string li2R = std::string("pCC#mu = ")  + VHS::ToString(pCCmu,7,5);
  std::string li3R = std::string("pCC#tau = ") + VHS::ToString(pCCtau,7,5);
  std::string li4R = std::string("pNC  = ")    + VHS::ToString(pNC,7,5);
  ptR->AddText(li1R.c_str());
  ptR->AddText(li2R.c_str());
  ptR->AddText(li3R.c_str());
  ptR->AddText(li4R.c_str());
  ptR->Draw();

  TPaveText* pt = new TPaveText(0.02,0.45,0.98,0.98,"brNDC");
  pt->SetLineWidth(2);
  pt->SetBorderSize(1);
  pt->SetTextSize(0.05);
  pt->SetTextAlign(12);

  std::string li1 =
    std::string("Event #") + VHS::ToString(nEntry+1) +
    std::string(" of ") + VHS::ToString(kNumEntries);
  std::string li2 = std::string("");
  if (evtt == vhsNC   ) li2 = "NC event ";
  if (evtt == vhsCCe  ) li2 = "e CCevent ";
  if (evtt == vhsCCmu ) li2 = "#mu CC event ";
  if (evtt == vhsCCtau) li2 = "#tau CC event ";
  li2 += (MC)?
    std::string("(") + VHS::ToString(NumPln) + std::string(" total planes)") :
    VHS::ToString(NumPln) + std::string(" total planes") ;
  std::string li3 = std::string("E_{vis} = ") + VHS::ToString(GeV,7,5) +
    std::string(" GeV (") + VHS::ToString(Sigcor,7,7) + std::string(" pe)");
  std::string li4 = (MC)?
    std::string("E_{#nu} = ") +
    VHS::ToString(MCenu,7,5)  +
    std::string(" GeV")       :
    "";
  std::string li5 = "";
  if ( dNC < dCCmu )
    li5 = std::string("Determined: NC event");
  if ( dCCmu < dNC )
    li5 = std::string("Determined: #mu CC event");
  //if ( dNC < dCCe && dNC < dCCmu && dNC < dCCtau )
  //li5 = std::string("Determined: NC event");
  //if ( dCCe < dNC && dCCe < dCCmu && dCCe < dCCtau )
  //li5 = std::string("Determined: e CC event");
  //if ( dCCmu < dNC && dCCmu < dCCe && dCCmu < dCCtau )
  //li5 = std::string("Determined: #mu CC event");
  //if ( dCCtau < dNC && dCCtau < dCCe && dCCtau < dCCmu )
  //li5 = std::string("Determined: #tau CC event");
  li5 += std::string(" (distance) ");
  pt->AddText(li1.c_str());
  pt->AddText(li2.c_str());
  pt->AddText(li3.c_str());
  if (MC) pt->AddText(li4.c_str());
  if ( (MC) ){
    if ((evtt == vhsNC    &&
         dNC > dCCmu       ) ||
        (evtt == vhsCCmu    &&
         dCCmu > dNC       )  )
         /*
         ( dNC    > dCCe  || dNC    > dCCmu || dNC    > dCCtau )) ||
        (evtt == vhsCCe   &&
         ( dCCe   > dNC   || dCCe   > dCCmu || dCCe   > dCCtau )) ||
        (evtt == vhsCCmu  &&
         ( dCCmu  > dNC   || dCCmu  > dCCe  || dCCmu  > dCCtau )) ||
        (evtt == vhsCCtau &&
         ( dCCtau > dNC   || dCCtau > dCCe  || dCCtau > dCCmu  ))  )
         */
      pt->AddText(li5.c_str())->SetTextColor(2);
    else
      pt->AddText(li5.c_str());
  } else {
    pt->AddText(li5.c_str());
  }

  std::string li6 = "";
  if ( pNC > pCCmu )
    li6 = std::string("Determined: NC event");
  if ( pCCmu > pNC )
    li6 = std::string("Determined: #mu CC event");
  /*
  if ( pNC > pCCe && pNC > pCCmu && pNC > pCCtau )
    li6 = std::string("Determined: NC event");
  if ( pCCe > pNC && pCCe > pCCmu && pCCe > pCCtau )
    li6 = std::string("Determined: e CC event");
  if ( pCCmu > pNC && pCCmu > pCCe && pCCmu > pCCtau )
    li6 = std::string("Determined: #mu CC event");
  if ( pCCtau > pNC && pCCtau > pCCe && pCCtau > pCCmu )
    li6 = std::string("Determined: #tau CC event");
  */
  li6 += std::string(" (loglikelihood) ");
  if ( (MC) ) {
    if ((evtt == vhsNC    &&
         pNC < pCCmu       ) ||
        (evtt == vhsCCmu  &&
         pCCmu < pNC       )  )
         /*
         ( pNC    < pCCe  || pNC    < pCCmu || pNC    < pCCtau )) ||
        (evtt == vhsCCe   &&
         ( pCCe   < pNC   || pCCe   < pCCmu || pCCe   < pCCtau )) ||
        (evtt == vhsCCmu  &&
         ( pCCmu  < pNC   || pCCmu  < pCCe  || pCCmu  < pCCtau )) ||
        (evtt == vhsCCtau &&
         ( pCCtau < pNC   || pCCtau < pCCe  || pCCtau < pCCmu  ))  )
         */
      pt->AddText(li6.c_str())->SetTextColor(2);
    else
      pt->AddText(li6.c_str());
  } else {
    pt->AddText(li6.c_str());
  }

  pt->Draw();

  delete evt; evt = 0;

  evtDisp->cd(3);

  evtDisp->SetEditable(kFALSE);

  return true;

}

//......................................................................

void VHS::FillDiscriminants(NtpSREvent*          evt,
                            TClonesArray*        stp,
                            std::vector<double>  avgNC,
                            std::vector<double>  avgCCe,
                            std::vector<double>  avgCCmu,
                            std::vector<double>  avgCCtau,
                            std::vector<double>  varNC,
                            std::vector<double>  varCCe,
                            std::vector<double>  varCCmu,
                            std::vector<double>  varCCtau,
                            std::vector<double>  pNC,
                            std::vector<double>  pCCe,
                            std::vector<double>  pCCmu,
                            std::vector<double>  pCCtau,
                            const int            nPlanes,
                            const int            nStrips,
                            const bool           bUnit,
                            VHSevent*&           vhsevt)
{

  // Gather our image and index sparse vectors
  std::vector<double> image;
  std::vector<int>    index;
  std::vector<double> theta;
  VHS::GetImage(evt->stp, evt->nstrip,
                stp, nPlanes, nStrips, image, index, theta);

  vhsevt->stp = image;
  vhsevt->ind = index;

  vhsevt->MOIThetaU = theta[0];
  vhsevt->MOIThetaV = theta[1];

  // Get minimum plane in U & V views
  vhsevt->MinPlnU = evt->plane.begu;
  vhsevt->MinPlnV = evt->plane.begv;

  // Fill relative loglikelihoods
  std::vector<double> fullImg = FullVec(image,index,nPlanes,nStrips);

  vhsevt->pNC    =
    VHS::GetLL(fullImg, pNC   , avgNC   , varNC   );
  vhsevt->pCCe   =
    VHS::GetLL(fullImg, pCCe  , avgCCe  , varCCe  );
  vhsevt->pCCmu  =
    VHS::GetLL(fullImg, pCCmu , avgCCmu , varCCmu );
  vhsevt->pCCtau =
    VHS::GetLL(fullImg, pCCtau, avgCCtau, varCCtau);

  // Fill mean distances
  if (bUnit) {
    VHS::UnitVector(fullImg );
    VHS::UnitVector(avgNC   );
    VHS::UnitVector(avgCCe  );
    VHS::UnitVector(avgCCmu );
    VHS::UnitVector(avgCCtau);
  }

  vhsevt->dNC    = VHS::GetDistance( fullImg, avgNC    );
  vhsevt->dCCe   = VHS::GetDistance( fullImg, avgCCe   );
  vhsevt->dCCmu  = VHS::GetDistance( fullImg, avgCCmu  );
  vhsevt->dCCtau = VHS::GetDistance( fullImg, avgCCtau );

}

//......................................................................

std::vector<double> VHS::FindMedian(std::vector< std::vector<double> > pt,
                                    std::vector< double >              avg,
                                    bool bUnit,
                                    bool bVerbose                         )
{
  const unsigned int nDim = avg.size(); // number of dimensions each point lives in
  const unsigned int nPts = pt.size();  // number of points for median-finding
  const double kEpsilon   = 1.e-6;      // precision to zero in calculations
  const double kDelta     = 1.e-4;      // cutoff for iterative median displacement

  if (bUnit) {
    VHS::UnitVector(avg);
    for (unsigned int i = 0; i < nPts; i++)
      VHS::UnitVector(pt[i]);
  }

  // Find the median via Weiszfeld iteration
  vector<double> median = avg;
  vector<double> lastMedian(nDim, 0.);
  double dDist    = 999.; // iterative difference in avg distance to iterative median
  double lastAvgD = 0.;   // average distance to the last iterative median
  double avgD     = 0.;   // average distance to the current iterative median

  // Find the average distance to our initial guess point (average)
  for (unsigned int i = 0; i < nPts; i++) {
    double dPt = VHS::GetDistance(median, pt[i]);
    avgD += dPt/nPts;
  }
  if (bVerbose)
    cout << "-- Iteration 0 : dItr = N/A : avgDist = " << avgD << " del(N/A)" << endl;

  // Iteratively move closer to the true median point
  double dItr     = 999.; // distance between last iterative median and new one
  int    nItr     = 0;    // number of iterations
  while( dItr > kDelta ) {
    lastMedian = median;
    for (unsigned int j = 0; j < nDim; j++) median[j] = 0.;
    lastAvgD = avgD;
    avgD = 0.;

    // Calculate our new median guess
    double norm = 0.;
    for (unsigned int i = 0; i < nPts; i++) {
      double dPt = VHS::GetDistance(lastMedian, pt[i]);

      if (!bUnit) norm += (dPt > kEpsilon)? 1./dPt : 1./kEpsilon;
      
      for (unsigned int j = 0; j < nDim; j++)
        median[j] += (dPt > kEpsilon)?
          pt[i][j]/dPt : pt[i][j]/TMath::Abs(kEpsilon) ;
    } // end for loop
    // Normalize our iterative median properly
    if (bUnit) VHS::UnitVector(median);
    else for (unsigned int j = 0; j < nDim; j++) median[j] *= 1./norm;

    dItr = VHS::GetDistance( lastMedian, median );
    // Find the average distance to our new guess point
    for (unsigned int i = 0; i < nPts; i++) {
      double dPt = VHS::GetDistance(median, pt[i]);
      avgD += dPt/nPts;
    }

    dDist = lastAvgD-avgD;
    nItr++;
    if (bVerbose)
      cout << "-- Iteration " << nItr << " : dItr = " << dItr
           << " : avgDist = " << avgD << " del(" << dDist
           << ")" << endl;
  } // end while
  if (bVerbose) {
    double dMedAvg = VHS::GetDistance( median, avg );
    cout << "|median-average| = " << dMedAvg << endl;
  }

  return median;

}

//......................................................................

std::vector<double> VHS::FullVec(std::vector<double> image,
                                 std::vector<int>    index,
                                 int                 nPlanes,
                                 int                 nStrips)
{
  std::vector<double> full( 2*nPlanes*nStrips, 0.);

  for (unsigned int i = 0; i < index.size(); i++)
    full[ index[i] ] = image[i];

  return full;

}

//......................................................................

double VHS::GetDistance(std::vector<double> vec0, std::vector<double> vec1)
{
  // Return the Euclidean distance between the points described
  // by the input vectors

  if (vec0.size() != vec1.size()) return -1.;

  double dist = 0.;

  for (unsigned int i = 0; i < vec0.size(); i++)
    dist += TMath::Power( vec0[i]-vec1[i], 2. );

  return TMath::Sqrt( dist );

}

//......................................................................

VHS::evtType VHS::GetEvtType( int inu, int iaction )
{
  if ( iaction == 0 ) return vhsNC; 

  if ( iaction == 1 && TMath::Abs(inu) == 12 ) return vhsCCe; 

  if ( iaction == 1 && TMath::Abs(inu) == 14 ) return vhsCCmu; 

  if ( iaction == 1 && TMath::Abs(inu) == 16 ) return vhsCCtau; 

  return vhsUnknown;

}

//......................................................................

void VHS::GetImage(int*                 index,
                   int                  nstp,
                   TClonesArray*        stp,
                   int                  nPlanes,
                   int                  nStrips,
                   std::vector<double>& image,
                   std::vector<int>&    vecInd,
                   std::vector<double>& theta)
{
  // Return a standard vector of hit values representing
  // the event image for both U & V views

  // The images will be centered on strip ~nStrips/2 & plane ~nPlanes/4

  image.clear();    // empty the return vector of values
  vecInd.clear();   // empty the return vector of indices
  theta.clear();    // empty the return vector of event angle
  theta = std::vector<double>(2,-999.);

  // Separate the input index vector & stp array into u & v
  std::vector< std::vector<int> > uvindex(2); //uvindex index 0=u, 1=v
  VHS::SeparateViews(index,nstp,stp,uvindex[0],uvindex[1]);

  const int kOffsetInd = nPlanes*nStrips;

  for (unsigned int uv = 0; uv <= 1 ; uv++) {
    std::vector<int> ind = uvindex[uv];

    // Derive image rotation info
    std::vector<int> centerPix(2,0);
    VHS::GetThetaAxis(stp, uvindex[uv], theta[uv], centerPix);

    // Find the minimum plane in the view
    int minPln    = 999;
    for (unsigned int j = 0; j < ind.size(); j++) {
      NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
      if (stpEntry) {
        int plane  = stpEntry->plane;
        int strip  = stpEntry->strip;
        VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
        if (plane < minPln) minPln = plane;
      } // end if stp ok
    } // end loop over strips

    // Find the end plane in our window
    int endPln = minPln + 2*nPlanes;

    // Find the weighted average plane in the view
    double stpXph    = 0.;
    double totSigcor = 0.;
    for (unsigned int j = 0; j < ind.size(); j++) {
      NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
      if (stpEntry) {
        int plane  = stpEntry->plane;
        int strip  = stpEntry->strip;
        VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
        if (plane >= minPln && plane < endPln) {
          double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;

          stpXph    += plane*sigcor;
          totSigcor += sigcor;
        } // end if strip is in our plane window
      } // end if stp ok
    } // end loop over strips

    // Find the weighted mean plane# in our window (~minos convention)
    // --> adjust the min & end plane numbers accordingly
    const int avgPln = TMath::FloorNint(stpXph/totSigcor);
    minPln = avgPln - TMath::FloorNint(0.25*2.*nPlanes);
    endPln = minPln + 2*nPlanes;

    // Find the weighted average strip in the view
    stpXph    = 0.;
    totSigcor = 0.;
    for (unsigned int j = 0; j < ind.size(); j++) {
      NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
      if (stpEntry) {
        int plane  = stpEntry->plane;
        int strip  = stpEntry->strip;
        VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
        if (plane >= minPln && plane < endPln) {
          double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;

          stpXph    += strip*sigcor;
          totSigcor += sigcor;
        } // end if strip is in our plane window
      } // end if stp ok
    } // end loop over strips

    // Find the weighted mean strip in our window
    const int avgStp = TMath::FloorNint(stpXph/totSigcor);
    const int minStp = avgStp - nStrips/2;

    // Loop over strips
    for (unsigned int i = 0; i < ind.size(); i++) {
      NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[i]));
      if (stpEntry) {
        int plane  = stpEntry->plane;
        int strip  = stpEntry->strip;
        VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
        double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;

        int iPln = (plane - minPln)/2;           // img coords
        int iStp = (strip - minStp);             // img coords
        int iVec = VecIndex(iPln,iStp,nPlanes);  // vec img coords

        if (iPln >= 0 && iPln < nPlanes &&
            iStp >= 0 && iStp < nStrips  ) {
          image.push_back(sigcor);
          if (uv == 0)
            vecInd.push_back(iVec);
          else
            vecInd.push_back(iVec+kOffsetInd);
        } // end if hit placement OK
      } // end if strip info OK
    } // end loop over strips
  } // end loop over u & v views

  return;

}

//......................................................................

double VHS::GetLL(std::vector<double> fullImage,
                  std::vector<double> pHit,
                  std::vector<double> avg,
                  std::vector<double> var)
{
  // Given a set of averages and variances,
  // calculate the Gaussian probability of the collection of
  // particular hit values convoluted with the probability of the
  // particular hit pattern

  double ll = 0.;
  const double kCutoff   = 1.e-10;
  const double kLnCutoff = TMath::Log10(kCutoff);

  for (unsigned int i = 0; i < fullImage.size(); i++) {
    double x     = fullImage[i];
    double mean  = avg[i];
    double sigma = TMath::Sqrt( var[i] );

    if ( x < kCutoff || pHit[i] < kCutoff ) 
      ll += TMath::Log10(1.-pHit[i]);
    else {
      ll += TMath::Log10(pHit[i]);
      ll += (TMath::Gaus(x,mean,sigma,true) > kCutoff)?
        TMath::Log10(TMath::Gaus(x,mean,sigma,true)) : kLnCutoff;
      ll -= (TMath::Gaus(mean,mean,sigma,true) > kCutoff)?
        TMath::Log10(TMath::Gaus(mean,mean,sigma,true)) : kLnCutoff;
    } // end if values within ranges
  } // end loop over image entries

  return ll;

}

//......................................................................

int VHS::GetPlane(int vecInd, int nPlanes)
{
  // Returns relative plane # (*not* MINOS global plane convention)
  return (nPlanes != 0)? vecInd % nPlanes : -1;
}

//......................................................................

int VHS::GetStrip(int vecInd, int nPlanes)
{
  // Returns relative strip # (*not* MINOS global strip convention)
  return (nPlanes != 0)? vecInd / nPlanes : -1;
}

//......................................................................

void VHS::GetThetaAxis(TClonesArray*     stp,
                       std::vector<int>  index,
                       double&           theta,
                       std::vector<int>& center)
{
  const unsigned int kNhits = index.size();

  // Find the weighted average plane & strip in the view
  // And the minimum plane/strip
  double stpXph    = 0.;
  double plnXph    = 0.;
  double totSigcor = 0.;
  int    minPln    = 999;
  int    minPlnStp = 999;
  for (unsigned int j = 0; j < kNhits; j++) {
    NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(index[j]));
    if (stpEntry) {
      int    plane  = stpEntry->plane;
      int    strip  = stpEntry->strip;
      double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;

      plnXph    += plane*sigcor;
      stpXph    += strip*sigcor;
      totSigcor += sigcor;

      if (plane < minPln) {
        minPln    = plane;
        minPlnStp = strip;
      } // end if minimum plane number
    } // end if stp ok
  } // end loop over strips

  center[0] = (int)TMath::Floor(plnXph/totSigcor);
  center[1] = (int)TMath::Floor(stpXph/totSigcor);

  double relPln = plnXph/totSigcor - 1.*minPln;
  double relStp = stpXph/totSigcor - 1.*minPlnStp;

  theta = (TMath::Abs(relPln) > 1.e-10)?
    TMath::ATan( relStp/relPln ) : 0.;

  return;
}

//......................................................................

std::vector<std::string> VHS::ParseNmList(const char* cstr)
{
  // Turn the input space delimited cstring into a vector of strings

  std::string strList = cstr;

  std::vector<std::string> strVec;

  int len = strList.length();

  unsigned int lastLoc = 0;
  unsigned int loc = strList.find(" ",lastLoc);
  while (lastLoc != std::string::npos) {
    int subLen = (loc != std::string::npos)? loc-lastLoc : len-lastLoc;
    std::string sub = strList.substr(lastLoc, subLen);
    if (sub.length() > 0) strVec.push_back(sub);
    lastLoc = (loc != std::string::npos)? loc+1 : loc;
    loc = strList.find(" ",lastLoc);
  }

  return strVec;
}

//......................................................................

void VHS::ReadFile(TFile*                inFile,
                   std::vector<double>&  avgNC,
                   std::vector<double>&  avgCCe,
                   std::vector<double>&  avgCCmu,
                   std::vector<double>&  avgCCtau,
                   std::vector<double>&  varNC,
                   std::vector<double>&  varCCe,
                   std::vector<double>&  varCCmu,
                   std::vector<double>&  varCCtau,
                   std::vector<double>&  pNChit,
                   std::vector<double>&  pCCehit,
                   std::vector<double>&  pCCmuhit,
                   std::vector<double>&  pCCtauhit,
                   const int             nPlanes)
{
  //
  // inFile must be open and valid
  //
  // Translate averages & completed variances from TH2D and write
  // to vectors for use
  //

  // Prepare to write our training info to output file
  inFile->cd();

  // NC hit prob hist
  const std::string UpNCName  = "UpNCHit";
  const std::string VpNCName  = "VpNCHit";
  TH2D* UpNCHist = (TH2D*)inFile->Get(UpNCName.c_str());
  TH2D* VpNCHist = (TH2D*)inFile->Get(VpNCName.c_str());
  std::vector<double> UpNChit = VHS::ToVector( UpNCHist, nPlanes);
  std::vector<double> VpNChit = VHS::ToVector( VpNCHist, nPlanes);
  pNChit.clear();
  for (unsigned int i = 0; i < UpNChit.size(); i++)
    pNChit.push_back(UpNChit[i]);
  for (unsigned int i = 0; i < VpNChit.size(); i++)
    pNChit.push_back(VpNChit[i]);

  // CCe hit prob hist
  const std::string UpCCeName  = "UpCCeHit";
  const std::string VpCCeName  = "VpCCeHit";
  TH2D* UpCCeHist = (TH2D*)inFile->Get(UpCCeName.c_str());
  TH2D* VpCCeHist = (TH2D*)inFile->Get(VpCCeName.c_str());
  std::vector<double> UpCCehit = VHS::ToVector( UpCCeHist, nPlanes);
  std::vector<double> VpCCehit = VHS::ToVector( VpCCeHist, nPlanes);
  pCCehit.clear();
  for (unsigned int i = 0; i < UpCCehit.size(); i++)
    pCCehit.push_back(UpCCehit[i]);
  for (unsigned int i = 0; i < VpCCehit.size(); i++)
    pCCehit.push_back(VpCCehit[i]);

  // CCmu hit prob hist
  const std::string UpCCmuName  = "UpCCmuHit";
  const std::string VpCCmuName  = "VpCCmuHit";
  TH2D* UpCCmuHist = (TH2D*)inFile->Get(UpCCmuName.c_str());
  TH2D* VpCCmuHist = (TH2D*)inFile->Get(VpCCmuName.c_str());
  std::vector<double> UpCCmuhit = VHS::ToVector( UpCCmuHist, nPlanes);
  std::vector<double> VpCCmuhit = VHS::ToVector( VpCCmuHist, nPlanes);
  pCCmuhit.clear();
  for (unsigned int i = 0; i < UpCCmuhit.size(); i++)
    pCCmuhit.push_back(UpCCmuhit[i]);
  for (unsigned int i = 0; i < VpCCmuhit.size(); i++)
    pCCmuhit.push_back(VpCCmuhit[i]);

  // CCtau hit prob hist
  const std::string UpCCtauName  = "UpCCtauHit";
  const std::string VpCCtauName  = "VpCCtauHit";
  TH2D* UpCCtauHist = (TH2D*)inFile->Get(UpCCtauName.c_str());
  TH2D* VpCCtauHist = (TH2D*)inFile->Get(VpCCtauName.c_str());
  std::vector<double> UpCCtauhit = VHS::ToVector( UpCCtauHist, nPlanes);
  std::vector<double> VpCCtauhit = VHS::ToVector( VpCCtauHist, nPlanes);
  pCCtauhit.clear();
  for (unsigned int i = 0; i < UpCCtauhit.size(); i++)
    pCCtauhit.push_back(UpCCtauhit[i]);
  for (unsigned int i = 0; i < VpCCtauhit.size(); i++)
    pCCtauhit.push_back(VpCCtauhit[i]);

  // NC Avg hist
  const std::string UavgNCName  = "UAvgNCImg";
  const std::string VavgNCName  = "VAvgNCImg";
  TH2D* UavgNCHist = (TH2D*)inFile->Get(UavgNCName.c_str());
  TH2D* VavgNCHist = (TH2D*)inFile->Get(VavgNCName.c_str());
  std::vector<double> UavgNC = VHS::ToVector( UavgNCHist, nPlanes);
  std::vector<double> VavgNC = VHS::ToVector( VavgNCHist, nPlanes);
  avgNC.clear();
  for (unsigned int i = 0; i < UavgNC.size(); i++)
    avgNC.push_back(UavgNC[i]);
  for (unsigned int i = 0; i < VavgNC.size(); i++)
    avgNC.push_back(VavgNC[i]);

  // CCe Avg hist
  const std::string UavgCCeName  = "UAvgCCeImg";
  const std::string VavgCCeName  = "VAvgCCeImg";
  TH2D* UavgCCeHist = (TH2D*)inFile->Get(UavgCCeName.c_str());
  TH2D* VavgCCeHist = (TH2D*)inFile->Get(VavgCCeName.c_str());
  std::vector<double> UavgCCe = VHS::ToVector( UavgCCeHist, nPlanes);
  std::vector<double> VavgCCe = VHS::ToVector( VavgCCeHist, nPlanes);
  avgCCe.clear();
  for (unsigned int i = 0; i < UavgCCe.size(); i++)
    avgCCe.push_back(UavgCCe[i]);
  for (unsigned int i = 0; i < VavgCCe.size(); i++)
    avgCCe.push_back(VavgCCe[i]);

  // CCmu Avg hist
  const std::string UavgCCmuName  = "UAvgCCmuImg";
  const std::string VavgCCmuName  = "VAvgCCmuImg";
  TH2D* UavgCCmuHist = (TH2D*)inFile->Get(UavgCCmuName.c_str());
  TH2D* VavgCCmuHist = (TH2D*)inFile->Get(VavgCCmuName.c_str());
  std::vector<double> UavgCCmu = VHS::ToVector( UavgCCmuHist, nPlanes);
  std::vector<double> VavgCCmu = VHS::ToVector( VavgCCmuHist, nPlanes);
  avgCCmu.clear();
  for (unsigned int i = 0; i < UavgCCmu.size(); i++)
    avgCCmu.push_back(UavgCCmu[i]);
  for (unsigned int i = 0; i < VavgCCmu.size(); i++)
    avgCCmu.push_back(VavgCCmu[i]);

  // CCtau Avg hist
  const std::string UavgCCtauName  = "UAvgCCtauImg";
  const std::string VavgCCtauName  = "VAvgCCtauImg";
  TH2D* UavgCCtauHist = (TH2D*)inFile->Get(UavgCCtauName.c_str());
  TH2D* VavgCCtauHist = (TH2D*)inFile->Get(VavgCCtauName.c_str());
  std::vector<double> UavgCCtau = VHS::ToVector( UavgCCtauHist, nPlanes);
  std::vector<double> VavgCCtau = VHS::ToVector( VavgCCtauHist, nPlanes);
  avgCCtau.clear();
  for (unsigned int i = 0; i < UavgCCtau.size(); i++)
    avgCCtau.push_back(UavgCCtau[i]);
  for (unsigned int i = 0; i < VavgCCtau.size(); i++)
    avgCCtau.push_back(VavgCCtau[i]);

  // NC variance hist
  const std::string UNCVarName  = "UNCVar";
  const std::string VNCVarName  = "VNCVar";
  TH2D* UNCVarHist = (TH2D*)inFile->Get(UNCVarName.c_str());
  TH2D* VNCVarHist = (TH2D*)inFile->Get(VNCVarName.c_str());
  std::vector<double> UNCVar = VHS::ToVector( UNCVarHist, nPlanes);
  std::vector<double> VNCVar = VHS::ToVector( VNCVarHist, nPlanes);
  varNC.clear();
  for (unsigned int i = 0; i < UNCVar.size(); i++)
    varNC.push_back(UNCVar[i]);
  for (unsigned int i = 0; i < VpNChit.size(); i++)
    varNC.push_back(VNCVar[i]);

  // CCe variance hist
  const std::string UCCeVarName  = "UCCeVar";
  const std::string VCCeVarName  = "VCCeVar";
  TH2D* UCCeVarHist = (TH2D*)inFile->Get(UCCeVarName.c_str());
  TH2D* VCCeVarHist = (TH2D*)inFile->Get(VCCeVarName.c_str());
  std::vector<double> UCCeVar = VHS::ToVector( UCCeVarHist, nPlanes);
  std::vector<double> VCCeVar = VHS::ToVector( VCCeVarHist, nPlanes);
  varCCe.clear();
  for (unsigned int i = 0; i < UCCeVar.size(); i++)
    varCCe.push_back(UCCeVar[i]);
  for (unsigned int i = 0; i < VpCCehit.size(); i++)
    varCCe.push_back(VCCeVar[i]);

  // CCmu variance hist
  const std::string UCCmuVarName  = "UCCmuVar";
  const std::string VCCmuVarName  = "VCCmuVar";
  TH2D* UCCmuVarHist = (TH2D*)inFile->Get(UCCmuVarName.c_str());
  TH2D* VCCmuVarHist = (TH2D*)inFile->Get(VCCmuVarName.c_str());
  std::vector<double> UCCmuVar = VHS::ToVector( UCCmuVarHist, nPlanes);
  std::vector<double> VCCmuVar = VHS::ToVector( VCCmuVarHist, nPlanes);
  varCCmu.clear();
  for (unsigned int i = 0; i < UCCmuVar.size(); i++)
    varCCmu.push_back(UCCmuVar[i]);
  for (unsigned int i = 0; i < VpCCmuhit.size(); i++)
    varCCmu.push_back(VCCmuVar[i]);

  // CCtau variance hist
  const std::string UCCtauVarName  = "UCCtauVar";
  const std::string VCCtauVarName  = "VCCtauVar";
  TH2D* UCCtauVarHist = (TH2D*)inFile->Get(UCCtauVarName.c_str());
  TH2D* VCCtauVarHist = (TH2D*)inFile->Get(VCCtauVarName.c_str());
  std::vector<double> UCCtauVar = VHS::ToVector( UCCtauVarHist, nPlanes);
  std::vector<double> VCCtauVar = VHS::ToVector( VCCtauVarHist, nPlanes);
  varCCtau.clear();
  for (unsigned int i = 0; i < UCCtauVar.size(); i++)
    varCCtau.push_back(UCCtauVar[i]);
  for (unsigned int i = 0; i < VpCCtauhit.size(); i++)
    varCCtau.push_back(VCCtauVar[i]);

  inFile->Close();

  return;

}

//......................................................................

void VHS::RotatePixel(int&         plane,    int&      strip,
                      const int    avgPlane, const int avgStrip,
                      const double theta                       )
{
  // rotate around pixel center to new coords (pln, stp);
  // --> round floats to nearest integer

  if ( TMath::Abs(theta) < 1.e-10 ) return; // ~zero rotation

  const double cosTheta = TMath::Cos(-theta); // rotate opposite dir as given
  const double sinTheta = TMath::Sin(-theta); // rotate opposite dir as given

  double newPlnD = cosTheta*(plane-avgPlane) + sinTheta*(strip-avgStrip);
  double newStpD = cosTheta*(strip-avgStrip) - sinTheta*(plane-avgPlane);

  if (newStpD-TMath::Floor(newStpD) < 0.5) newStpD = TMath::Floor(newStpD);
  else newStpD = TMath::Ceil(newStpD);

  if (newPlnD-TMath::Floor(newPlnD) < 0.5) newPlnD = TMath::Floor(newPlnD);
  else newPlnD = TMath::Ceil(newPlnD);

  strip = avgStrip + (int)newStpD;
  plane = avgPlane + (int)newPlnD;

  return;

}

//......................................................................

void VHS::SeparateViews(int* index, int nstp, TClonesArray* stp,
                        std::vector<int>& ustp, std::vector<int>& vstp)
{
  // Sort hit strips into separate views

  for (int i = 0; i < nstp; i++) {
    int ind = index[i];
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>(stp->At(ind));
    if (strip) {
      unsigned int pv = strip->planeview;
      if (pv == PlaneView::kU) ustp.push_back(ind);
      if (pv == PlaneView::kV) vstp.push_back(ind);
    }
  }
  return;
}

//......................................................................

std::vector< VHSevent* > VHS::Skim(NtpStRecord*         ntpst,
                                   std::vector<double>  avgNC,
                                   std::vector<double>  avgCCe,
                                   std::vector<double>  avgCCmu,
                                   std::vector<double>  avgCCtau,
                                   std::vector<double>  varNC,
                                   std::vector<double>  varCCe,
                                   std::vector<double>  varCCmu,
                                   std::vector<double>  varCCtau,
                                   std::vector<double>  pNC,
                                   std::vector<double>  pCCe,
                                   std::vector<double>  pCCmu,
                                   std::vector<double>  pCCtau,
                                   const int            nPlanes,
                                   const int            nStrips,
                                   const bool           bUnit)
{
  // Inputs: NtpStRecord, average images, variances,
  //         vector of raw hit probabilities, window size (pln x stp),
  //         flag signalling whether to compute distances on unit or
  //         absolute image vectors
  // Output: Vector of pointers to complete VHSevent records
  //

  // Define the vector we will return
  vector< VHSevent* > vhslist;

  // If we have no record, do nothing
  if ( !ntpst ) return vhslist;

  const NtpSREventSummary& evthdr = ntpst->evthdr;
  TClonesArray*            evt    = ntpst->evt;
  TClonesArray*            stp    = ntpst->stp;
  TClonesArray*            mc     = ntpst->mc;
  TClonesArray*            thevt  = ntpst->thevt;

  // Event quantities
  if (evt && stp) {
    const unsigned int nEvt = evthdr.nevent; // evt->GetEntries();

    for (unsigned int i = 0; i < nEvt; i++){
      NtpSREvent* evtEntry = dynamic_cast<NtpSREvent*>( evt->At(i) );
      if ( !evtEntry ) continue;

      // Define our VHSevent to fill
      VHSevent* vhsevt = new VHSevent( ntpst->GetHeader() ) ;
      vhsevt->SetName("VHSevent");
      vhsevt->SetTitle("VHS Event");

      // Fill NtpSREvent info
      vhsevt->GeV    = evtEntry->ph.gev;
      vhsevt->NumPln = evtEntry->plane.end - evtEntry->plane.beg;
      vhsevt->NumStp = evtEntry->nstrip;
      vhsevt->Sigcor = evtEntry->ph.sigcor;

      VHS::evtType evtt = vhsUnknown;

      // Fill MC information
      if ( mc && thevt ) {
        NtpTHEvent*  thevtEntry  = dynamic_cast<NtpTHEvent*>(thevt->At(i));
        NtpMCTruth*  mcEntry     = (thevtEntry)?
          dynamic_cast<NtpMCTruth*>( mc->At( (thevtEntry->neumc) ) ) : 0;
        if (mcEntry) {
          vhsevt->MCenu     = mcEntry->p4neu[3];
          vhsevt->MCinu     = mcEntry->inu;
          vhsevt->MCiaction = mcEntry->iaction;
          vhsevt->MCrescode = mcEntry->iresonance;
          vhsevt->MCQsq     = mcEntry->q2;
          vhsevt->MCWsq     = mcEntry->w2;
          vhsevt->MCx       = mcEntry->x;
          vhsevt->MCy       = mcEntry->y;

          evtt = VHS::GetEvtType(mcEntry->inu, mcEntry->iaction);
        } // end if mcEntry
      } // end if mc && thevt
      vhsevt->MC = (evtt != vhsUnknown);

      VHS::FillDiscriminants(evtEntry, stp,
                             avgNC, avgCCe, avgCCmu, avgCCtau,
                             varNC, varCCe, varCCmu, varCCtau,
                             pNC,   pCCe,   pCCmu,   pCCtau,
                             nPlanes, nStrips, bUnit, vhsevt );

      vhslist.push_back( vhsevt ); // add this event to the list

    } // end loop over event entries
  } // end if evt and stp arrays exist

  return vhslist;
}

//......................................................................

std::vector<double> VHS::SubtractVec(std::vector<double> image0,
                                     std::vector<double> image1)
{
  // 'nuff said

  std::vector<double> subt(image0.size(), 0.);

  if (image0.size() != image1.size() ) return subt;

  for (unsigned int i = 0; i < image0.size(); i++)
    subt[i] = image0[i] - image1[i];

  return subt;

}

//......................................................................

template<class T>
std::string VHS::ToString(const T& thing, int w, int p)
{

  std::ostringstream os;
  os << std::setw(w) << std::setprecision(p) << thing;
  return os.str();

}

//......................................................................

void VHS::ToTH2D(std::vector<double> vec,
                 int         nPlanes, int         nStrips,
                 std::string name,    std::string title,
                 TH2D*&      Uview,   TH2D*&      Vview  )
{
  // 'nuff said

  std::string Uname  = "U" + name;
  std::string Utitle = "U" + title;

  std::string Vname  = "V" + name;
  std::string Vtitle = "V" + title;

  //if (Uview) delete Uview;
  //if (Vview) delete Vview;

  Uview = new TH2D(Uname.c_str(), Utitle.c_str(),
                   nPlanes, 0, nPlanes,
                   nStrips, 0, nStrips);
  Vview = new TH2D(Vname.c_str(), Vtitle.c_str(),
                   nPlanes, 0, nPlanes,
                   nStrips, 0, nStrips);
  const unsigned int kOffset = nPlanes*nStrips;
  for (unsigned int i = 0; i < vec.size(); i++){
    if ( i < kOffset )
      Uview->SetBinContent( VHS::GetPlane(i,nPlanes)+1,
                            VHS::GetStrip(i,nPlanes)+1,
                            vec[i]                  );
    else
      Vview->SetBinContent( VHS::GetPlane(i-kOffset,nPlanes)+1,
                            VHS::GetStrip(i-kOffset,nPlanes)+1,
                            vec[i]                  );
  } // end loop over vec entries

  return;

}

//......................................................................

std::vector<double> VHS::ToVector(TH2D* hist, int nPlanes)
{
  // 'nuff said

  const int nPln  = hist->GetNbinsX();
  const int nStp  = hist->GetNbinsY();
  const int nBins = nPln*nStp;

  std::vector<double> vec(nBins,0.);

  if (nPln != nPlanes) return vec;

  for (int ipln = 0; ipln < nPln; ipln++)
    for (int istp = 0; istp < nStp; istp++)
      vec[VecIndex(ipln,istp,nPlanes)] =
        hist->GetBinContent(ipln+1, istp+1);

  return vec;

}

//......................................................................

int  VHS::Train(NtpStRecord*          ntpst,
                std::vector<double>&  avgNC,
                std::vector<double>&  avgCCe,
                std::vector<double>&  avgCCmu,
                std::vector<double>&  avgCCtau,
                std::vector<double>&  varNC,
                std::vector<double>&  varCCe,
                std::vector<double>&  varCCmu,
                std::vector<double>&  varCCtau,
                std::vector<int>&     numNC,
                std::vector<int>&     numCCe,
                std::vector<int>&     numCCmu,
                std::vector<int>&     numCCtau,
                int&                  NCevts,
                int&                  eCCevts,
                int&                  muCCevts,
                int&                  tauCCevts,
                const int             maxTrain,
                const int             nPlanes,
                const int             nStrips,
                const int             cutPlanes,
                const double          eRecoMax,
                const double          eTrueMax)
{
  // Input : NtpStRecord, average images to modify, variances to
  //         modify, vectors of number of hits to each pixel to
  //         modify, #planes in each view window, #strips in each view
  //         window
  // Output: Modified average images & modified covariance matrices
  //

  int nEvtUsed = 0;

  // If we have no record, do nothing
  if ( !ntpst ) return nEvtUsed;

  const NtpSREventSummary& evthdr = ntpst->evthdr;
  TClonesArray*            evt    = ntpst->evt;
  TClonesArray*            stp    = ntpst->stp;
  TClonesArray*            mc     = ntpst->mc;
  TClonesArray*            thevt  = ntpst->thevt;

  // Event quantities
  if (evt && stp && mc && thevt) {
    const unsigned int nEvt = evthdr.nevent; // evt->GetEntries();

    for (unsigned int i = 0; i < nEvt; i++){
      NtpSREvent* evtEntry = dynamic_cast<NtpSREvent*>( evt->At(i) );
      if ( !evtEntry ) continue;

      NtpTHEvent*  thevtEntry  = dynamic_cast<NtpTHEvent*>(thevt->At(i));
      NtpMCTruth*  mcEntry     = (thevtEntry)?
        dynamic_cast<NtpMCTruth*>( mc->At( (thevtEntry->neumc) ) ) : 0;
      if ( !mcEntry ) continue;

      const double eReco = evtEntry->ph.gev;
      const double eTrue = mcEntry->p4neu[3];

      if (eReco < eRecoMax || eTrue < eTrueMax ) continue;

      VHS::evtType evtt = VHS::GetEvtType(mcEntry->inu, mcEntry->iaction);

      if ( evtt == vhsUnknown ) continue;

      if (evtEntry->plane.end - evtEntry->plane.beg > cutPlanes) continue;

      if ( ( evtt == vhsNC    && NCevts    >= maxTrain ) ||
           ( evtt == vhsCCe   && eCCevts   >= maxTrain ) ||
           ( evtt == vhsCCmu  && muCCevts  >= maxTrain ) ||
           ( evtt == vhsCCtau && tauCCevts >= maxTrain )  )
        continue;

      // Gather our image and index sparse vectors
      std::vector<double> image;
      std::vector<int>    index;
      std::vector<double> theta;
      VHS::GetImage(evtEntry->stp, evtEntry->nstrip,
                    stp, nPlanes, nStrips, image, index, theta);

      // Break up into cases by event type
      switch (evtt) {
        case vhsNC :
          for (unsigned int ind = 0; ind < image.size(); ind++) {
            int    ii = index[ind];
            double xi = image[ind];
            int    ni = numNC[ii];
            
            avgNC[ii] = (avgNC[ii]*ni + xi   )/(ni+1);
            varNC[ii] = (varNC[ii]*ni + xi*xi)/(ni+1);
            numNC[ii]++;
          } // end outer for loop 
          NCevts++;

          break;
        case vhsCCe :
          for (unsigned int ind = 0; ind < image.size(); ind++) {
            int    ii = index[ind];
            double xi = image[ind];
            int    ni = numCCe[ii];
            
            avgCCe[ii] = (avgCCe[ii]*ni + xi   )/(ni+1);
            varCCe[ii] = (varCCe[ii]*ni + xi*xi)/(ni + 1);
            numCCe[ii]++;
          } // end outer for loop 
          eCCevts++;

          break;
        case vhsCCmu :
          for (unsigned int ind = 0; ind < image.size(); ind++) {
            int    ii = index[ind];
            double xi = image[ind];
            int    ni = numCCmu[ii];
            
            avgCCmu[ii] = (avgCCmu[ii]*ni + xi   )/(ni+1);
            varCCmu[ii] = (varCCmu[ii]*ni + xi*xi)/(ni+1);
            numCCmu[ii]++;
          } // end outer for loop 
          muCCevts++;

          break;
        case vhsCCtau :
          for (unsigned int ind = 0; ind < image.size(); ind++) {
            int    ii = index[ind];
            double xi = image[ind];
            int    ni = numCCtau[ii];
            
            avgCCtau[ii] = (avgCCtau[ii]*ni + xi   )/(ni+1);
            varCCtau[ii] = (varCCtau[ii]*ni + xi*xi)/(ni+1);
            numCCtau[ii]++;
          } // end outer for loop 
          tauCCevts++;

          break;
        case vhsUnknown :
          continue;
      } // end switch (evtType)

      nEvtUsed++;

    } // end for loop over events
  } // end if the ntpst arrays are valid

  return nEvtUsed;

}

//......................................................................

void VHS::TrainPost(TFile*                outFile,
                    std::vector<double>&  avgNC,
                    std::vector<double>&  avgCCe,
                    std::vector<double>&  avgCCmu,
                    std::vector<double>&  avgCCtau,
                    std::vector<double>&  varNC,
                    std::vector<double>&  varCCe,
                    std::vector<double>&  varCCmu,
                    std::vector<double>&  varCCtau,
                    std::vector<int>&     numNC,
                    std::vector<int>&     numCCe,
                    std::vector<int>&     numCCmu,
                    std::vector<int>&     numCCtau,
                    int                   NCevts,
                    int                   eCCevts,
                    int                   muCCevts,
                    int                   tauCCevts,
                    const int             nPlanes,
                    const int             nStrips)
{
  //
  // Post-processing of the variances, assuming input "variances"
  // are entries simply of <x^2>
  //
  // Translate averages & completed variances into TH2D and write
  // to file, outFile
  //
  // Call this only once at the end of the total training job
  //

  const unsigned int kImgSz   = avgNC.size();

  // Calculate the simple raw hit probabilities
  std::vector<double> pNChit   (kImgSz, 0.0);
  std::vector<double> pCCehit  (kImgSz, 0.0);
  std::vector<double> pCCmuhit (kImgSz, 0.0);
  std::vector<double> pCCtauhit(kImgSz, 0.0);

  for (unsigned int pix = 0; pix < kImgSz; pix++) {
    pNChit[pix]    = (NCevts > 0)?
      (double)numNC[pix]    / (double)NCevts    : 0.;
    pCCehit[pix]   = (eCCevts > 0)?
      (double)numCCe[pix]   / (double)eCCevts   : 0.;
    pCCmuhit[pix]  = (muCCevts > 0)?
      (double)numCCmu[pix]  / (double)muCCevts  : 0.;
    pCCtauhit[pix] = (tauCCevts > 0)?
      (double)numCCtau[pix] / (double)tauCCevts : 0.;
  }

  // Complete the variance calculation: var = <x^2> - <x>^2
  for (unsigned int ii = 0; ii < kImgSz; ii++) {
    varNC[ii]    -= avgNC[ii]    * avgNC[ii];
    varCCe[ii]   -= avgCCe[ii]   * avgCCe[ii];
    varCCmu[ii]  -= avgCCmu[ii]  * avgCCmu[ii];
    varCCtau[ii] -= avgCCtau[ii] * avgCCtau[ii];

  } // end strip loop

  VHS::WriteFile( outFile,
                  avgNC, avgCCe, avgCCmu, avgCCtau,
                  varNC, varCCe, varCCmu, varCCtau,
                  pNChit, pCCehit, pCCmuhit, pCCtauhit,
                  nPlanes, nStrips);

  return;

}

//......................................................................

void VHS::UnitVector(std::vector<double>& vec)
{
  double len = 0.;

  for (unsigned int i = 0; i < vec.size(); i++)
    len += TMath::Power( vec[i], 2. );

  len = TMath::Sqrt( len );

  if (len <= 0.) return;

  for (unsigned int i = 0; i < vec.size(); i++)
    vec[i] = vec[i] / len ;

  return;

}

//......................................................................

int VHS::VecIndex(int ipln, int istp, int nPlanes)
{
  // A particular choice of convention

  return (ipln + nPlanes * istp);

}

//......................................................................

void VHS::WriteFile(TFile*                outFile,
                    std::vector<double>&  avgNC,
                    std::vector<double>&  avgCCe,
                    std::vector<double>&  avgCCmu,
                    std::vector<double>&  avgCCtau,
                    std::vector<double>&  varNC,
                    std::vector<double>&  varCCe,
                    std::vector<double>&  varCCmu,
                    std::vector<double>&  varCCtau,
                    std::vector<double>&  pNChit,
                    std::vector<double>&  pCCehit,
                    std::vector<double>&  pCCmuhit,
                    std::vector<double>&  pCCtauhit,
                    const int             nPlanes,
                    const int             nStrips)
{
  //
  // outFile must be open and valid
  //
  // Translate averages & completed variances into TH2D and write
  // to file, outFile
  //
  // Call this only once at the end of the total training job
  //

  // Prepare to write our training info to output file
  outFile->cd();

  const std::string imgXTitle = "Plane";
  const std::string imgYTitle = "Strip";

  // NC hit prob hist
  const std::string pNCTitle = "NC Hit Probabilities";
  const std::string pNCName  = "pNCHit";
  TH2D* UpNCHist; TH2D* VpNCHist;
  VHS::ToTH2D( pNChit,
               nPlanes, nStrips,
               pNCName.c_str(),
               pNCTitle.c_str(),
               UpNCHist, VpNCHist);
  UpNCHist->SetXTitle( imgXTitle.c_str() );
  UpNCHist->SetYTitle( imgYTitle.c_str() );
  UpNCHist->Write();
  VpNCHist->SetXTitle( imgXTitle.c_str() );
  VpNCHist->SetYTitle( imgYTitle.c_str() );
  VpNCHist->Write();

  // CCe hit prob hist
  const std::string pCCeTitle = "CCe Hit Probabilities";
  const std::string pCCeName  = "pCCeHit";
  TH2D* UpCCeHist; TH2D* VpCCeHist;
  VHS::ToTH2D( pCCehit,
               nPlanes, nStrips,
               pCCeName.c_str(),
               pCCeTitle.c_str(),
               UpCCeHist, VpCCeHist);
  UpCCeHist->SetXTitle( imgXTitle.c_str() );
  UpCCeHist->SetYTitle( imgYTitle.c_str() );
  UpCCeHist->Write();
  VpCCeHist->SetXTitle( imgXTitle.c_str() );
  VpCCeHist->SetYTitle( imgYTitle.c_str() );
  VpCCeHist->Write();

  // CCmu hit prob hist
  const std::string pCCmuTitle = "CC#mu Hit Probabilities";
  const std::string pCCmuName  = "pCCmuHit";
  TH2D* UpCCmuHist; TH2D* VpCCmuHist;
  VHS::ToTH2D( pCCmuhit,
               nPlanes, nStrips,
               pCCmuName.c_str(),
               pCCmuTitle.c_str(),
               UpCCmuHist, VpCCmuHist);
  UpCCmuHist->SetXTitle( imgXTitle.c_str() );
  UpCCmuHist->SetYTitle( imgYTitle.c_str() );
  UpCCmuHist->Write();
  VpCCmuHist->SetXTitle( imgXTitle.c_str() );
  VpCCmuHist->SetYTitle( imgYTitle.c_str() );
  VpCCmuHist->Write();

  // CCtau hit prob hist
  const std::string pCCtauTitle = "CC#tau Hit Probabilities";
  const std::string pCCtauName  = "pCCtauHit";
  TH2D* UpCCtauHist; TH2D* VpCCtauHist;
  VHS::ToTH2D( pCCtauhit,
               nPlanes, nStrips,
               pCCtauName.c_str(),
               pCCtauTitle.c_str(),
               UpCCtauHist, VpCCtauHist);
  UpCCtauHist->SetXTitle( imgXTitle.c_str() );
  UpCCtauHist->SetYTitle( imgYTitle.c_str() );
  UpCCtauHist->Write();
  VpCCtauHist->SetXTitle( imgXTitle.c_str() );
  VpCCtauHist->SetYTitle( imgYTitle.c_str() );
  VpCCtauHist->Write();

  // NC Avg hist
  const std::string avgNCTitle = "Average NC Image";
  const std::string avgNCName  = "AvgNCImg";
  TH2D* UavgNCHist; TH2D* VavgNCHist;
  VHS::ToTH2D( avgNC,
               nPlanes, nStrips,
               avgNCName.c_str(),
               avgNCTitle.c_str(),
               UavgNCHist, VavgNCHist);
  UavgNCHist->SetXTitle( imgXTitle.c_str() );
  UavgNCHist->SetYTitle( imgYTitle.c_str() );
  UavgNCHist->Write();
  VavgNCHist->SetXTitle( imgXTitle.c_str() );
  VavgNCHist->SetYTitle( imgYTitle.c_str() );
  VavgNCHist->Write();

  // CCe Avg hist
  const std::string avgCCeTitle = "Average CCe Image";
  const std::string avgCCeName  = "AvgCCeImg";
  TH2D* UavgCCeHist; TH2D* VavgCCeHist;
  VHS::ToTH2D( avgCCe,
               nPlanes, nStrips,
               avgCCeName.c_str(),
               avgCCeTitle.c_str(),
               UavgCCeHist, VavgCCeHist);
  UavgCCeHist->SetXTitle( imgXTitle.c_str() );
  UavgCCeHist->SetYTitle( imgYTitle.c_str() );
  UavgCCeHist->Write();
  VavgCCeHist->SetXTitle( imgXTitle.c_str() );
  VavgCCeHist->SetYTitle( imgYTitle.c_str() );
  VavgCCeHist->Write();

  // CCmu Avg hist
  const std::string avgCCmuTitle = "Average CC#mu Image";
  const std::string avgCCmuName  = "AvgCCmuImg";
  TH2D* UavgCCmuHist; TH2D* VavgCCmuHist;
  VHS::ToTH2D( avgCCmu,
               nPlanes, nStrips,
               avgCCmuName.c_str(),
               avgCCmuTitle.c_str(),
               UavgCCmuHist, VavgCCmuHist);
  UavgCCmuHist->SetXTitle( imgXTitle.c_str() );
  UavgCCmuHist->SetYTitle( imgYTitle.c_str() );
  UavgCCmuHist->Write();
  VavgCCmuHist->SetXTitle( imgXTitle.c_str() );
  VavgCCmuHist->SetYTitle( imgYTitle.c_str() );
  VavgCCmuHist->Write();

  // CCtau Avg hist
  const std::string avgCCtauTitle = "Average CC#tau Image";
  const std::string avgCCtauName  = "AvgCCtauImg";
  TH2D* UavgCCtauHist; TH2D* VavgCCtauHist;
  VHS::ToTH2D( avgCCtau,
               nPlanes, nStrips,
               avgCCtauName.c_str(),
               avgCCtauTitle.c_str(),
               UavgCCtauHist, VavgCCtauHist);
  UavgCCtauHist->SetXTitle( imgXTitle.c_str() );
  UavgCCtauHist->SetYTitle( imgYTitle.c_str() );
  UavgCCtauHist->Write();
  VavgCCtauHist->SetXTitle( imgXTitle.c_str() );
  VavgCCtauHist->SetYTitle( imgYTitle.c_str() );
  VavgCCtauHist->Write();

  // NC variance hist
  const std::string varNCTitle = "NC Variance";
  const std::string varNCName  = "NCVar";
  TH2D* UvarNCHist; TH2D* VvarNCHist;
  VHS::ToTH2D( varNC,
               nPlanes, nStrips,
               varNCName.c_str(),
               varNCTitle.c_str(),
               UvarNCHist, VvarNCHist);
  UvarNCHist->SetXTitle( imgXTitle.c_str() );
  UvarNCHist->SetYTitle( imgYTitle.c_str() );
  UvarNCHist->Write();
  VvarNCHist->SetXTitle( imgXTitle.c_str() );
  VvarNCHist->SetYTitle( imgYTitle.c_str() );
  VvarNCHist->Write();

  // CCe variance hist
  const std::string varCCeTitle = "CCe Variance";
  const std::string varCCeName  = "CCeVar";
  TH2D* UvarCCeHist; TH2D* VvarCCeHist;
  VHS::ToTH2D( varCCe,
               nPlanes, nStrips,
               varCCeName.c_str(),
               varCCeTitle.c_str(),
               UvarCCeHist, VvarCCeHist);
  UvarCCeHist->SetXTitle( imgXTitle.c_str() );
  UvarCCeHist->SetYTitle( imgYTitle.c_str() );
  UvarCCeHist->Write();
  VvarCCeHist->SetXTitle( imgXTitle.c_str() );
  VvarCCeHist->SetYTitle( imgYTitle.c_str() );
  VvarCCeHist->Write();

  // CCmu variance hist
  const std::string varCCmuTitle = "CCmu Variance";
  const std::string varCCmuName  = "CCmuVar";
  TH2D* UvarCCmuHist; TH2D* VvarCCmuHist;
  VHS::ToTH2D( varCCmu,
               nPlanes, nStrips,
               varCCmuName.c_str(),
               varCCmuTitle.c_str(),
               UvarCCmuHist, VvarCCmuHist);
  UvarCCmuHist->SetXTitle( imgXTitle.c_str() );
  UvarCCmuHist->SetYTitle( imgYTitle.c_str() );
  UvarCCmuHist->Write();
  VvarCCmuHist->SetXTitle( imgXTitle.c_str() );
  VvarCCmuHist->SetYTitle( imgYTitle.c_str() );
  VvarCCmuHist->Write();

  // CCtau variance hist
  const std::string varCCtauTitle = "CCtau Variance";
  const std::string varCCtauName  = "CCtauVar";
  TH2D* UvarCCtauHist; TH2D* VvarCCtauHist;
  VHS::ToTH2D( varCCtau,
               nPlanes, nStrips,
               varCCtauName.c_str(),
               varCCtauTitle.c_str(),
               UvarCCtauHist, VvarCCtauHist);
  UvarCCtauHist->SetXTitle( imgXTitle.c_str() );
  UvarCCtauHist->SetYTitle( imgYTitle.c_str() );
  UvarCCtauHist->Write();
  VvarCCtauHist->SetXTitle( imgXTitle.c_str() );
  VvarCCtauHist->SetYTitle( imgYTitle.c_str() );
  VvarCCtauHist->Write();

  return;

}

//......................................................................

---