NuExtraction.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Jul/2005       
//                                                                    
// Contact: j.j.hartnell@rl.ac.uk

#include <set>
#include <cmath>

#include "TClonesArray.h"
#include "TH2F.h"
#include "TMath.h"
#include "TVector3.h"

#include "BeamDataUtil/BDSpillAccessor.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "DcsUser/CoilTools.h"
#include "DataUtil/DataQualDB.h"
#include "DataUtil/EnergyCorrections.h"
#include "MessageService/MsgFormat.h"
#include "CandFitTrackSA/Ntp/NtpFitSARecord.h"
#include "CandFitTrackSA/Ntp/NtpFitSAPlane.h"
#include "CandFitTrackSA/Ntp/NtpFitSA.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRTimeStatus.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "StandardNtuple/NtpStRecord.h"

#include "SpillTiming/SpillTimeFinder.h"
#include "SpillTiming/SpillServerMonFinder.h"
#include "RunQuality/DbuFarRunQuality.h"
#include "RunQuality/DbuNearRunQuality.h"

#include "NtupleUtils/LISieve.h"
#include "NtupleUtils/NuLibrary.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuExtraction.h"
#include "NtupleUtils/NuMCEvent.h"

// for the NC analysis
#include "AnalysisNtuples/Module/CondensedNtpModule.h"
#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "AnalysisNtuples/Module/ANtpInfoObjectFillerNC.h"
#include "AnalysisNtuples/ANtpEventInfoNC.h"
#include "AnalysisNtuples/ANtpTrackInfoNC.h"
#include "AnalysisNtuples/ANtpShowerInfoNC.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AstroUtil/AstTime.h"

using std::endl;
using std::cout;
using std::map;
using std::vector;

namespace EnergyCorrections {}
using namespace EnergyCorrections;

CVSID("$Id: NuExtraction.cxx,v 1.67 2010/01/25 13:17:43 evansj Exp $");

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

NuExtraction::NuExtraction()
{
  MSG("NuExtraction",Msg::kDebug)
    <<"Running NuExtraction Constructor..."<<endl;


  MSG("NuExtraction",Msg::kDebug)
    <<"Finished NuExtraction Constructor"<<endl;
}

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

NuExtraction::~NuExtraction()
{
  MSG("NuExtraction",Msg::kDebug)
    <<"Running NuExtraction Destructor..."<<endl;
  

  MSG("NuExtraction",Msg::kDebug)
    <<"Finished NuExtraction Destructor"<<endl;
}

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

void NuExtraction::ExtractBasicInfo(const NtpStRecord& ntp,
                                    const NtpSREvent& evt,
                                    NuEvent& nu) const
{

  //get the run, snarl, etc info
  this->ExtractGeneralInfo(ntp,nu);
  //get info from the evt
  this->ExtractEvtInfo(evt,nu);
  //extract trk info (needed to get best track info)
  this->ExtractTrkInfo(ntp,evt,nu);
  //extract shw info
  this->ExtractShwInfo(ntp,evt,nu);     
}

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

void NuExtraction::ExtractAuxiliaryInfo(const NtpStRecord& ntp,
                                        const NtpSREvent& evt,
                                        const NtpFitSARecord* pntpSA,
                                        NuEvent& nu) const
{
  //extract a bunch of supplementary information needed for analysis
  
  //get the SA fitter information if it is available
  this->ExtractSAFitInfo(pntpSA,nu);
  //get the timing info
  this->ExtractMinMaxEvtTimes(ntp,evt,nu);
  this->ExtractTimeToNearestSpill(nu);
  //get the data quality
  this->ExtractDataQuality(nu);
  this->ExtractLITags(ntp,nu);
  this->ExtractBeamInfoDB(nu);//done per snarl but do here as safer
  this->ExtractAdditionalDSTVariables(nu);
}

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

void NuExtraction::ExtractGeneralInfo(const NtpStRecord& ntp,
                                      NuEvent& nu) const
{
  const RecCandHeader& rec=ntp.GetHeader();

  nu.run=rec.GetRun();
  nu.subRun=rec.GetSubRun();
  nu.runType=rec.GetRunType();
  nu.errorCode=rec.GetErrorCode();
  nu.snarl=rec.GetSnarl();
  nu.trigSrc=rec.GetTrigSrc();
  nu.timeFrame=rec.GetTimeFrame();
  nu.remoteSpillType=rec.GetRemoteSpillType();

  nu.detector=rec.GetVldContext().GetDetector();
  nu.simFlag=rec.GetVldContext().GetSimFlag();
  nu.timeSec=rec.GetVldContext().GetTimeStamp().GetSec();
  nu.timeNanoSec=rec.GetVldContext().GetTimeStamp().GetNanoSec();
  nu.timeSeconds=rec.GetVldContext().GetTimeStamp().GetSeconds();
  
  const NtpSREventSummary& evthdr=ntp.evthdr;
  nu.trigtime=evthdr.trigtime;
  nu.planeEvtHdrBeg=evthdr.plane.beg;
  nu.planeEvtHdrEnd=evthdr.plane.end;
  nu.snarlPulseHeight=evthdr.ph.sigcor;

  //a hack to correct the trigger bit of the cosmic MC in the ND
  if (nu.simFlag==SimFlag::kMC && nu.detector==Detector::kNear) {
    //the cosmic MC doesn't have any entries
    TClonesArray& mcTca=*ntp.mc;
    if (mcTca.GetEntries()==0){ 
      MAXMSG("MeuAnalysis",Msg::kWarning,5)
        <<"Performing hack to correct the fTrigSrc of"
        <<" ND cosmic MC, was="
        <<nu.trigSrc<<", setting to 4"<<endl;
      nu.trigSrc=4;//a plane trigger
    }
  }
}

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

void NuExtraction::ExtractSAFitInfo(const NtpFitSARecord* pntpSA,
                                    NuEvent& nu) const
{
  if (!pntpSA) {
    MAXMSG("MeuAnalysis",Msg::kDebug,5)
      <<"ExtractSAFitInfo: no NtpFitSARecord"<<endl;
    return;
  }

  //get a reference
  const NtpFitSARecord& ntpSA=*pntpSA;
  const TClonesArray& fitsaTca=(*ntpSA.fitsa);

  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();
  
  if (nu.trkIndex1>=0) {
    const NtpFitSA& fitsa=
      *(dynamic_cast<NtpFitSA*>(fitsaTca[nu.trkIndex1]));
    
    nu.trkfitpassSA1=fitsa.fit.pass;
    nu.trkvtxdcoszSA1=fitsa.fit.dcosz; 
    nu.chargeSA1=lib.reco.GetTrackCharge(fitsa.fit.qp, nu.hornIsReverse);
    nu.qpSA1=fitsa.fit.qp;
    nu.sigqpSA1=fitsa.fit.eqp;
    nu.chi2SA1=fitsa.fit.chi2;
    nu.ndofSA1=fitsa.fit.ndf;
    nu.probSA1=TMath::Prob(fitsa.fit.chi2,fitsa.fit.ndf);
    TVector3 xyz=lib.reco.uvz2xyz(fitsa.fit.u,fitsa.fit.v,fitsa.fit.z,
                                  lib.reco.GetVldContext(nu));
    nu.xTrkVtxSA1=xyz.X();
    nu.yTrkVtxSA1=xyz.Y();
    nu.zTrkVtxSA1=fitsa.fit.z;
    nu.uTrkVtxSA1=fitsa.fit.u;
    nu.vTrkVtxSA1=fitsa.fit.v;
  }

  if (nu.trkIndex2>=0) {
    const NtpFitSA& fitsa=
      *(dynamic_cast<NtpFitSA*>(fitsaTca[nu.trkIndex2]));
    
    nu.trkfitpassSA2=fitsa.fit.pass;
    nu.trkvtxdcoszSA2=fitsa.fit.dcosz; 
    nu.chargeSA2=lib.reco.GetTrackCharge(fitsa.fit.qp, nu.hornIsReverse);
    nu.qpSA2=fitsa.fit.qp;
    nu.sigqpSA2=fitsa.fit.eqp;
    nu.chi2SA2=fitsa.fit.chi2;
    nu.ndofSA2=fitsa.fit.ndf;
    nu.probSA2=TMath::Prob(fitsa.fit.chi2,fitsa.fit.ndf);
    TVector3 xyz=lib.reco.uvz2xyz(fitsa.fit.u,fitsa.fit.v,fitsa.fit.z,
                                  lib.reco.GetVldContext(nu));
    nu.xTrkVtxSA2=xyz.X();
    nu.yTrkVtxSA2=xyz.Y();
    nu.zTrkVtxSA2=fitsa.fit.z;
    nu.uTrkVtxSA2=fitsa.fit.u;
    nu.vTrkVtxSA2=fitsa.fit.v;
  }

  if (nu.trkIndex3>=0) {
    const NtpFitSA& fitsa=
      *(dynamic_cast<NtpFitSA*>(fitsaTca[nu.trkIndex3]));
    
    nu.trkfitpassSA3=fitsa.fit.pass;
    nu.trkvtxdcoszSA3=fitsa.fit.dcosz; 
    nu.chargeSA3=lib.reco.GetTrackCharge(fitsa.fit.qp, nu.hornIsReverse);
    nu.qpSA3=fitsa.fit.qp;
    nu.sigqpSA3=fitsa.fit.eqp;
    nu.chi2SA3=fitsa.fit.chi2;
    nu.ndofSA3=fitsa.fit.ndf;
    nu.probSA3=TMath::Prob(fitsa.fit.chi2,fitsa.fit.ndf);
    TVector3 xyz=lib.reco.uvz2xyz(fitsa.fit.u,fitsa.fit.v,fitsa.fit.z,
                                  lib.reco.GetVldContext(nu));
    nu.xTrkVtxSA3=xyz.X();
    nu.yTrkVtxSA3=xyz.Y();
    nu.zTrkVtxSA3=fitsa.fit.z;
    nu.uTrkVtxSA3=fitsa.fit.u;
    nu.vTrkVtxSA3=fitsa.fit.v;
  }

  //copy across the variables for the best track
  lib.reco.SetBestTrkSAFit(nu);
}
/*
  MAXMSG("MeuAnalysis",Msg::kInfo,500)
    <<"nu.ntrk="<<nu.ntrk
    <<" ntrk="<<ntpSA.ntrack
    <<", nu.planeTrkBeg="<<nu.planeTrkBeg
    <<", nu.planeTrkEnd="<<nu.planeTrkEnd
    <<endl;

  MAXMSG("MeuAnalysis",Msg::kInfo,500)
    <<"nu.trkIndex="<<nu.trkIndex<<endl;

  for (Int_t i=0;i<ntpSA.ntrack;i++) {      
    const NtpFitSA& fitsa=
      *(dynamic_cast<NtpFitSA*>(fitsaTca[i]));
    
    MAXMSG("MeuAnalysis",Msg::kInfo,500)
      <<"i="<<i
      <<", pln.begin="<<fitsa.plane.begin
      <<", pln.end="<<fitsa.plane.end
      <<"SA: pass="<<fitsa.fit.pass
      <<", qp="<<fitsa.fit.qp
      <<endl;
  }
  
  if (nu.trkIndex<0) {
    MAXMSG("MeuAnalysis",Msg::kWarning,50)
      <<"ExtractSAFitInfo: nu.trkIndex="<<nu.trkIndex<<endl;
  }


    MAXMSG("MeuAnalysis",Msg::kInfo,500)
      <<"SR: pass="<<nu.trkfitpass
      <<", qp="<<nu.qp
      <<", eqp="<<nu.sigqp
      <<", chi2="<<nu.chi2
      <<", ndof="<<nu.ndof
      <<", prob="<<nu.prob
      <<endl
      <<"SA: pass="<<fitsa.fit.pass
      <<", qp="<<fitsa.fit.qp
      <<", eqp="<<fitsa.fit.eqp
      <<", chi2="<<fitsa.fit.chi2
      <<", ndof="<<fitsa.fit.ndf
      <<", prob="<<TMath::Prob(fitsa.fit.chi2,fitsa.fit.ndf)
      <<endl;
      //chargeSA
  */

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

void NuExtraction::ExtractCoilInfo(NuEvent& nu) const
{
  //return if not data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting Coil info for simFlag="<<nu.simFlag<<endl;
    //set these to be the most likely values
    nu.coilIsOk=true;
    int fieldNo = NuUtilities::GetDigit(nu.run, 5);
    if (fieldNo == 1 || fieldNo == 3) {
      MAXMSG("NuExtraction",Msg::kInfo,1)
      << "Found field id #" << fieldNo << ", setting coil current to forward."
      <<endl;            
      nu.coilIsReverse=false;
    }
    else if (fieldNo == 2 || fieldNo == 4) {
      MAXMSG("NuExtraction",Msg::kInfo,1)
      << "Found field id #" << fieldNo << ", setting coil current to reverse."
      <<endl;                  
      nu.coilIsReverse=true;
    }
    else {
      MAXMSG("NuExtraction",Msg::kWarning,1)
      << "Don't recognize field id #" << fieldNo << ", only 1, 2, 3, 4 are allowed values.  Assuming field is forward."
      <<endl;      
      nu.coilIsReverse=false;
    }
    return;
  }
  
  //check whether to use the database
  if (!nu.useDBForDataQuality) {//use data quality flag
    MAXMSG("NuExtraction",Msg::kWarning,1)
      <<"Not extracting CoilCurrent (flag set to not query database)"
      <<endl;
    //leave values as current (for re-reconstruction case)
    return;
  }

  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();
  VldContext vc=lib.reco.GetVldContext(nu);
  
  //extract the information about the coil using CoilTools
  //this is much faster than using the database
  nu.coilIsOk=CoilTools::IsOK(vc);
  nu.coilIsReverse=CoilTools::IsReverse(vc);
  
  //get the coil current here too
  this->ExtractCoilCurrent(nu);
}

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

void NuExtraction::ExtractCoilCurrent(NuEvent& nu) const
{
  //return if not data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting CoilCurrent for simFlag="<<nu.simFlag<<endl;
    return;
  }
  
  //check whether to use the database
  if (!nu.useDBForDataQuality) {//use data quality flag
    MAXMSG("NuExtraction",Msg::kWarning,1)
      <<"Not extracting CoilCurrent (flag set to not query database)"
      <<endl;
    //leave values as they are (for re-reconstruction case)
    return;
  }

  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //get vc
  VldContext vc=lib.reco.GetVldContext(nu);

  //get the current, depending on detector
  if (nu.detector==Detector::kNear) {
    const Dcs_Mag_Near* magnear=CoilTools::Instance().GetMagNear(vc);
    if (magnear) nu.coilCurrent=magnear->GetCurrent();
    else {
      MAXMSG("NuExtraction",Msg::kWarning,10)
        <<"No rows in Dcs_Mag_Near DB table for VldContext="<<vc
        <<" setting coilCurrent=0"<<endl;
      nu.coilCurrent=0;//set to default value
    }
  } 
  else if (nu.detector==Detector::kFar) {
    const BfldDbiCoilState* coilstate1=
      CoilTools::Instance().GetCoilState(vc,1);
    const BfldDbiCoilState* coilstate2=
      CoilTools::Instance().GetCoilState(vc,2);
    //ahh, two values, one slot ... take average
    if (coilstate1 && coilstate2) {
      nu.coilCurrent=0.5*(coilstate1->GetCurrent()+
                          coilstate2->GetCurrent());
    } 
    else {
      MAXMSG("NuExtraction",Msg::kInfo,10)
        <<"Insufficient rows in BfldDbiCoilState DB table, VldContext="
        <<vc<<" setting coilCurrent=0"<<endl;
      nu.coilCurrent=0;//set to default value
    }
  }
  else cout<<"Ahhh, det="<<nu.detector<<endl;
  
  //The old slow way:
  //DbiResultPtr<Dcs_Mag_Near> ptr(vc);
  //if(ptr.GetNumRows()){
  //const Dcs_Mag_Near* row=ptr.GetRow(0);
  //nu.coilCurrent=row->GetCurrent();
  //}
  //else {
  //MAXMSG("NuExtraction",Msg::kInfo,10)
  //  <<"No rows in Dcs_Mag_Near DB table for VldContext="<<vc
  //  <<" setting coilCurrent=0"<<endl;
  //nu.coilCurrent=0;//set to default value
  //}

  MAXMSG("NuExtraction",Msg::kInfo,10)
    <<"Extracted CoilCurrent="<<nu.coilCurrent<<endl;
}

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

void NuExtraction::ExtractLITags(const NtpStRecord& ntp,NuEvent& nu) const
{
  
  //get an event header
  const NtpSREventSummary& evthdr=ntp.evthdr;      
  
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //extract all the possible LI variables
  nu.litime=evthdr.litime;
  nu.isLI=LISieve::IsLI(ntp);
  nu.litag=lib.reco.FDRCBoundary(nu);
}

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

void NuExtraction::ExtractBeamInfoDB(NuEvent& nu) const
{
  //check if data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting BeamInfoDB for simFlag="<<nu.simFlag
      <<endl;
    return;
  }

  //check whether to use the database
  if (!nu.useDBForBeamInfo) {
    MAXMSG("NuExtraction",Msg::kWarning,1)
      <<"Not extracting BeamInfoDB (flag set to not query database)"
      <<endl;
    //leave values as current (for re-reconstruction case)
    return;
  }

  MAXMSG("NuExtraction",Msg::kInfo,1)
    <<"Extracting BeamInfoDB for simFlag="<<nu.simFlag
    <<", detector="<<nu.detector<<endl;

  //do it per snarl (especially for the ND)
  MAXMSG("NuExtraction",Msg::kWarning,1)
    <<"ExtractBeamInfoDB: this needs to be sorted out properly"<<endl;

  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();
  
  //get the event time stamp
  VldContext vc=lib.reco.GetVldContext(nu);
  VldTimeStamp vts=vc.GetTimeStamp();
  
  //get a spill accessor
  BDSpillAccessor &sa = BDSpillAccessor::Get();
  const BeamMonSpill *bs = sa.LoadSpill(vts);

  //now determine if good beam
  Bool_t goodBeam=false;
  BMSpillAna bmana;
  bmana.UseDatabaseCuts();
  bmana.SetSpill(*bs);

  //set the time difference to determine if spill is in time
  VldTimeStamp delta_vts=vts-bs->SpillTime();
  bmana.SetTimeDiff(delta_vts.GetSeconds());
  
  //determine whether to select
  if(bmana.SelectSpill()) goodBeam=true;
  else goodBeam=false;
  
  //store the variable
  nu.goodBeam=goodBeam;
  
  // Extract run period for data from VldContext
  static VldTimeStamp* vtsRunIEnd = new VldTimeStamp(2006,4,1,0,0,0);
  static VldTimeStamp* vtsRunIIEnd = new VldTimeStamp(2007,8,1,0,0,0);
  static VldTimeStamp* vtsRunIIIEnd = new VldTimeStamp(2009,8,1,0,0,0);
  static VldTimeStamp* beginRHC = new VldTimeStamp(2009,9,29,0,0,0);
  static VldTimeStamp* endRHC = new VldTimeStamp(2010,3,1,0,0,0);
  
  MAXMSG("NuExtraction",Msg::kInfo, 1) << "beginRHC=" << beginRHC->AsString() << ", endRHC=" << endRHC->AsString() << endl;
  
  //store the beam variables (these are junk if the spill is not in time)
  nu.hornCur=bs->fHornCur;
  nu.beamTypeDB=bs->BeamType();
  nu.potDB=bs->fTortgt;

  if(vts.GetDate()<vtsRunIEnd->GetDate()) nu.runPeriod = 1;
  else if(vts.GetDate()<vtsRunIIEnd->GetDate()) nu.runPeriod = 2;
  else if(vts.GetDate()<vtsRunIIIEnd->GetDate()) nu.runPeriod = 3;
  else nu.runPeriod = 4;
  
  // Check for RHC
  if (vts.GetDate() >= beginRHC->GetDate() && 
      vts.GetDate() < endRHC->GetDate()  ) {
    MAXMSG("NuExtraction",Msg::kInfo, 10) << "Found RHC date: " << vts.AsString() << endl;
    nu.hornIsReverse =  true;
    nu.hornCur *= -1;
    if (nu.beamTypeDB == BeamType::kL010z185i) {
      nu.beamTypeDB = BeamType::kL010z185i_rev;
    }
    else {
      MAXMSG("NuExtraction",Msg::kInfo, 50) << "RHC event has non-rev beamtype: " 
      << BeamType::AsString(static_cast<BeamType::BeamType_t>(nu.beamTypeDB)) 
      << " (" << nu.beamTypeDB << ")" << endl;
    }
  }
  else {
    MAXMSG("NuExtraction",Msg::kInfo, 10) << "Found FHC date: " << vts.AsString() << endl;
    nu.hornIsReverse = false;
  }
  
  
  MAXMSG("NuExtraction",Msg::kInfo,10)
  <<"Extracting BeamInfoDB found" 
  <<  " goodBeam=" << nu.goodBeam
  << ", hornCur=" << nu.hornCur
  << ", beamTypeDB=" << nu.beamTypeDB
  << ", potDB=" << nu.potDB
  << ", hornIsReverse="<< (nu.hornIsReverse?"Yes":"No")
  << ", runPeriod=" << nu.runPeriod
  << endl;
  
  
  
  //other possibilities
  //Float_t hbw=bs->fProfWidX*1000.;//make it in mm
  //Float_t vbw=bs->fProfWidY*1000.;//make it in mm
  //Float_t hpos1=bs->fTargProfX*1000.;//make it in mm
  //Float_t vpos1=bs->fTargProfY*1000.;//make it in mm
  //BeamType::BeamType_t beamType=bs->BeamType();
  //Bool_t horn_on = bs->GetStatusBits().horn_on;
  //Bool_t target_in = bs->GetStatusBits().target_in;
}

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

void NuExtraction::ExtractTimeToNearestSpill(NuEvent& nu) const
{
  //check if data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting TimeToNearestSpill for simFlag="<<nu.simFlag
      <<endl;
    //set a default value
    nu.timeToNearestSpill=-999999;
    return;
  }

  //check if the far detector
  if (nu.detector!=Detector::kFar) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting TimeToNearestSpill for detector="<<nu.detector
      <<endl;
    //set a default value
    nu.timeToNearestSpill=-999999;
    return;
  }

  //check whether to use the database
  if (!nu.useDBForSpillTiming) {
    MAXMSG("NuExtraction",Msg::kWarning,1)
      <<"Not extracting TimeToNearestSpill"
      <<" (flag set to not query database)"
      <<", nu.timeToNearestSpill="<<nu.timeToNearestSpill
      <<endl;
    //leave as current value (for re-reconstruction case)
    return;
  }
 
  MAXMSG("NuExtraction",Msg::kInfo,1)
    <<"Extracting TimeToNearestSpill for simFlag="<<nu.simFlag
    <<", detector="<<nu.detector<<endl;
  SpillTimeFinder& spillFinder=SpillTimeFinder::Instance();
  
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  VldContext vc=lib.reco.GetVldContext(nu);
  nu.timeToNearestSpill=spillFinder.GetTimeToNearestSpill(vc);

  VldTimeStamp nearestspill = spillFinder.GetTimeOfNearestSpill(vc);
  nu.nearestSpillSec = nearestspill.GetSec();
  nu.nearestSpillNanosec = nearestspill.GetNanoSec();
}

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

void NuExtraction::ExtractDataQuality(NuEvent& nu) const
{  
  //return if not data
  if (nu.simFlag!=SimFlag::kData) {
    MAXMSG("NuExtraction",Msg::kInfo,1)
      <<"Not extracting Data Quality info for simFlag="<<nu.simFlag<<endl;
    //set these to be the most likely values
    nu.isGoodDataQuality = 1;
    nu.isGoodDataQualityRUN = 1;
    nu.isGoodDataQualityCOIL = 1;
    nu.isGoodDataQualityHV = 1;
    nu.isGoodDataQualityGPS = 1;
    return;
  }

  //check whether to use the database
  if (!nu.useDBForDataQuality) {
    MAXMSG("NuExtraction",Msg::kWarning,1)
      <<"Not extracting DataQuality"
      <<" (flag set to not query database)"
      <<", isGoodDataQuality="<<nu.isGoodDataQuality
      <<endl;
    //nu.isGoodDataQuality=leave as current value
    return;
  }
    
  // get validity context from ntuple record
  //const RecCandHeader* ntpHeader = &(ntp.GetHeader());
  //VldContext cx = ntpHeader->GetVldContext();
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();
  VldContext cx = lib.reco.GetVldContext(nu);
  
  // check data quality
  // ==================
  // run DataUtil methods

  nu.isGoodDataQuality = DataUtil::IsGoodData(cx);
  nu.isGoodDataQualityRUN = DataUtil::IsGoodDataRUN(cx);
  nu.isGoodDataQualityCOIL = DataUtil::IsGoodDataCOIL(cx);
  nu.isGoodDataQualityHV = DataUtil::IsGoodDataHV(cx);
  nu.isGoodDataQualityGPS = DataUtil::IsGoodDataGPS(cx);

  // extract run quality variables
  // =============================

  // far detector - search DbuFarRunQuality table
  if( cx.GetDetector() == Detector::kFar ) {   
    // perform database query
    DbiResultPtr<DbuFarRunQuality> farquery("DBUFARRUNQUALITY",cx,0);

    if( farquery.GetNumRows()>0 ){
      const DbuFarRunQuality* rowptr = farquery.GetRow(0);
      nu.numActiveCrates = rowptr->GetCrateMask();
      nu.numTimeFrames = rowptr->GetTimeFrames();
      nu.numGoodSnarls = rowptr->GetGoodSnarls();
      nu.snarlRateMedian = rowptr->GetMedianSnarlRate();
      nu.snarlRateMax = rowptr->GetMaxSnarlRate();
    }
  }

  // near detector
  if( cx.GetDetector() == Detector::kNear ) {
    //perform database query
    DbiResultPtr<DbuNearRunQuality> nearquery("DBUNEARRUNQUALITY",cx,0);
    
    if( nearquery.GetNumRows()>0 ){
      const DbuNearRunQuality* rowptr = nearquery.GetRow(0);
      nu.numActiveCrates = 8 - rowptr->GetColdCrates();
      nu.numTimeFrames = rowptr->GetSubrunLength();
      nu.numGoodSnarls = rowptr->GetTriggersSpill();
      nu.snarlRateMedian = rowptr->GetSnarlRateMedian();
      nu.snarlRateMax = rowptr->GetSnarlRateMax();
    }
  }


  // extract gps data
  // ================
  // use SpillServerMonFinder to access nearest spill monitoring block,
  // only necessary for far detector
  if( cx.GetDetector() == Detector::kFar ) {
  
    SpillServerMonFinder& spillFinder = SpillServerMonFinder::Instance();
    const SpillServerMon& nearestSpill = spillFinder.GetNearestSpill(cx);
    VldTimeStamp dt = nearestSpill.GetSpillTime()-cx.GetTimeStamp();

    nu.deltaSecToSpillGPS = dt.GetSec();
    nu.deltaNanoSecToSpillGPS = dt.GetNanoSec();
    nu.gpsError = nearestSpill.GetSpillTimeError();
    nu.gpsSpillType = nearestSpill.GetSpillType();
  }  

  // coil currents
  // =============
  // set this information in ExtractCoilInfo() method
}

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

void NuExtraction::ExtractTrkInfo(const NtpStRecord& ntp,
                                  const NtpSREvent& evt,
                                  NuEvent& nu) const
{
  Msg::LogLevel_t logLevel=Msg::kDebug;  
  
  MAXMSG("NuExtraction",logLevel,200)
    <<"evt="<<evt.index<<", trks="<<evt.ntrack
    <<", shws="<<evt.nshower<<endl;
  
  //extract the trk info
  this->ExtractPrimaryTrkInfo(ntp,evt,nu);
  this->ExtractSecondTrkInfo(ntp,evt,nu);
  this->ExtractThirdTrkInfo(ntp,evt,nu);
}

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

void NuExtraction::ExtractTrkShwInfo(const NtpStRecord& ntp,
                               const NtpSREvent& evt,
                               NuEvent& nu) const
{
  Msg::LogLevel_t logLevel=Msg::kDebug;  
  
  MAXMSG("NuExtraction",logLevel,200)
    <<"evt="<<evt.index<<", trks="<<evt.ntrack
    <<", shws="<<evt.nshower<<endl;
  
  //extract the trk info
  this->ExtractTrkInfo(ntp,evt,nu);
  
  //extract the shower info (does 1-3)
  this->ExtractShwInfo(ntp,evt,nu);
}

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

void NuExtraction::ExtractEvtInfo(const NtpSREvent& evt,
                            NuEvent& nu) const
{
    
  nu.evt=evt.index;
  nu.slc=evt.slc;
  nu.ndigitEvt=evt.ndigit;
  nu.nstripEvt=evt.nstrip;
  nu.nshw=evt.nshower;
  nu.ntrk=evt.ntrack;
  nu.primshw=evt.primshw;
  nu.primtrk=evt.primtrk;
  nu.rawPhEvt=evt.ph.raw;
  nu.evtphsigcor=evt.ph.sigcor;
  nu.evtphsigmap=evt.ph.sigmap;
  nu.planeEvtN=evt.plane.n;
  nu.planeEvtNu=evt.plane.nu;
  nu.planeEvtNv=evt.plane.nv;

  nu.xEvtVtx=evt.vtx.x;
  nu.yEvtVtx=evt.vtx.y;
  nu.zEvtVtx=evt.vtx.z;
  nu.uEvtVtx=evt.vtx.u;
  nu.vEvtVtx=evt.vtx.v;
  nu.planeEvtVtx=evt.vtx.plane;
  nu.planeEvtBeg=evt.plane.beg;
  nu.planeEvtBegu=evt.plane.begu;
  nu.planeEvtBegv=evt.plane.begv;
  
  nu.xEvtEnd=evt.end.x;
  nu.yEvtEnd=evt.end.y;
  nu.zEvtEnd=evt.end.z;
  nu.uEvtEnd=evt.end.u;
  nu.vEvtEnd=evt.end.v;
  nu.planeEvtEnd=evt.plane.end;
  nu.planeEvtEndu=evt.plane.endu;
  nu.planeEvtEndv=evt.plane.endv;
}

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

void NuExtraction::ExtractShwInfo(const NtpStRecord& ntp,
                            const NtpSREvent& evt,
                            NuEvent& nu) const
{
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //get pointer to shower
  const NtpSRShower* pshw=lib.reco.GetShowerWithIndexX(ntp,evt,0);
  nu.shwExists1=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists1=true;
    nu.ndigitShw1=shw.ndigit;
    nu.nstripShw1=shw.nstrip;
    nu.nplaneShw1=shw.plane.n;
    nu.shwEnCor1=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor1=shw.shwph.linCCgev;//has to go first
    nu.shwEnLinCCCor1=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnLinCCNoCor1=shw.shwph.linCCgev;
    nu.shwEnLinNCCor1=lib.reco.GetShowerEnergyCor(shw.shwph.linNCgev,CandShowerHandle::kNC,nu);
    nu.shwEnLinNCNoCor1=shw.shwph.linNCgev;
    nu.shwEnWtCCCor1=lib.reco.GetShowerEnergyCor(shw.shwph.wtCCgev,CandShowerHandle::kWtCC,nu);
    nu.shwEnWtCCNoCor1=shw.shwph.wtCCgev;
    nu.shwEnWtNCCor1=lib.reco.GetShowerEnergyCor(shw.shwph.wtNCgev,CandShowerHandle::kWtNC,nu);
    nu.shwEnWtNCNoCor1=shw.shwph.wtNCgev;

    nu.shwEnMip1=shw.ph.mip;//Note: shw.ph.gev=shw.shwph.wtCCgev
    nu.planeShwBeg1=shw.plane.beg;
    nu.planeShwEnd1=shw.plane.end;
    nu.planeShwMax1=lib.reco.GetShwMaxPlane(ntp, shw);
    nu.xShwVtx1=shw.vtx.x;
    nu.yShwVtx1=shw.vtx.y;
    nu.zShwVtx1=shw.vtx.z;
  }
  
  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,1);
  nu.shwExists2=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists2=true;
    nu.ndigitShw2=shw.ndigit;
    nu.nstripShw2=shw.nstrip;
    nu.nplaneShw2=shw.plane.n;
    nu.shwEnCor2=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor2=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnLinCCCor2=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnLinCCNoCor2=shw.shwph.linCCgev;
    nu.shwEnLinNCCor2=lib.reco.GetShowerEnergyCor(shw.shwph.linNCgev,CandShowerHandle::kNC,nu);
    nu.shwEnLinNCNoCor2=shw.shwph.linNCgev;
    nu.shwEnWtCCCor2=lib.reco.GetShowerEnergyCor(shw.shwph.wtCCgev,CandShowerHandle::kWtCC,nu);
    nu.shwEnWtCCNoCor2=shw.shwph.wtCCgev;
    nu.shwEnWtNCCor2=lib.reco.GetShowerEnergyCor(shw.shwph.wtNCgev,CandShowerHandle::kWtNC,nu);
    nu.shwEnWtNCNoCor2=shw.shwph.wtNCgev;
    nu.shwEnMip2=shw.ph.mip;
    nu.planeShwBeg2=shw.plane.beg;
    nu.planeShwEnd2=shw.plane.end;
    nu.planeShwMax2=lib.reco.GetShwMaxPlane(ntp, shw);
    nu.xShwVtx2=shw.vtx.x;
    nu.yShwVtx2=shw.vtx.y;
    nu.zShwVtx2=shw.vtx.z;
  }

  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,2);
  nu.shwExists3=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists3=true;
    nu.ndigitShw3=shw.ndigit;
    nu.nstripShw3=shw.nstrip;
    nu.nplaneShw3=shw.plane.n;
    nu.shwEnCor3=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor3=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnLinCCCor3=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnLinCCNoCor3=shw.shwph.linCCgev;
    nu.shwEnLinNCCor3=lib.reco.GetShowerEnergyCor(shw.shwph.linNCgev,CandShowerHandle::kNC,nu);
    nu.shwEnLinNCNoCor3=shw.shwph.linNCgev;
    nu.shwEnWtCCCor3=lib.reco.GetShowerEnergyCor(shw.shwph.wtCCgev,CandShowerHandle::kWtCC,nu);
    nu.shwEnWtCCNoCor3=shw.shwph.wtCCgev;
    nu.shwEnWtNCCor3=lib.reco.GetShowerEnergyCor(shw.shwph.wtNCgev,CandShowerHandle::kWtNC,nu);
    nu.shwEnWtNCNoCor3=shw.shwph.wtNCgev;
    nu.shwEnMip3=shw.ph.mip;
    nu.planeShwBeg3=shw.plane.beg;
    nu.planeShwEnd3=shw.plane.end;
    nu.planeShwMax3=lib.reco.GetShwMaxPlane(ntp, shw);
    nu.xShwVtx3=shw.vtx.x;
    nu.yShwVtx3=shw.vtx.y;
    nu.zShwVtx3=shw.vtx.z;
  }


  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,3);
  nu.shwExists4=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists4=true;
    nu.ndigitShw4=shw.ndigit;
    nu.nstripShw4=shw.nstrip;
    nu.nplaneShw4=shw.plane.n;
    nu.shwEnCor4=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor4=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnLinCCCor4=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnLinCCNoCor4=shw.shwph.linCCgev;
    nu.shwEnLinNCCor4=lib.reco.GetShowerEnergyCor(shw.shwph.linNCgev,CandShowerHandle::kNC,nu);
    nu.shwEnLinNCNoCor4=shw.shwph.linNCgev;
    nu.shwEnWtCCCor4=lib.reco.GetShowerEnergyCor(shw.shwph.wtCCgev,CandShowerHandle::kWtCC,nu);
    nu.shwEnWtCCNoCor4=shw.shwph.wtCCgev;
    nu.shwEnWtNCCor4=lib.reco.GetShowerEnergyCor(shw.shwph.wtNCgev,CandShowerHandle::kWtNC,nu);
    nu.shwEnWtNCNoCor4=shw.shwph.wtNCgev;
    nu.shwEnMip4=shw.ph.mip;
    nu.planeShwBeg4=shw.plane.beg;
    nu.planeShwEnd4=shw.plane.end;
    nu.planeShwMax4=lib.reco.GetShwMaxPlane(ntp, shw);
    nu.xShwVtx4=shw.vtx.x;
    nu.yShwVtx4=shw.vtx.y;
    nu.zShwVtx4=shw.vtx.z;
  }

  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,4);
  nu.shwExists5=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists5=true;
    nu.ndigitShw5=shw.ndigit;
    nu.nstripShw5=shw.nstrip;
    nu.nplaneShw5=shw.plane.n;
    nu.shwEnCor5=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor5=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnLinCCCor5=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnLinCCNoCor5=shw.shwph.linCCgev;
    nu.shwEnLinNCCor5=lib.reco.GetShowerEnergyCor(shw.shwph.linNCgev,CandShowerHandle::kNC,nu);
    nu.shwEnLinNCNoCor5=shw.shwph.linNCgev;
    nu.shwEnWtCCCor5=lib.reco.GetShowerEnergyCor(shw.shwph.wtCCgev,CandShowerHandle::kWtCC,nu);
    nu.shwEnWtCCNoCor5=shw.shwph.wtCCgev;
    nu.shwEnWtNCCor5=lib.reco.GetShowerEnergyCor(shw.shwph.wtNCgev,CandShowerHandle::kWtNC,nu);
    nu.shwEnWtNCNoCor5=shw.shwph.wtNCgev;
    nu.shwEnMip5=shw.ph.mip;
    nu.planeShwBeg5=shw.plane.beg;
    nu.planeShwEnd5=shw.plane.end;
    nu.planeShwMax5=lib.reco.GetShwMaxPlane(ntp, shw);
    nu.xShwVtx5=shw.vtx.x;
    nu.yShwVtx5=shw.vtx.y;
    nu.zShwVtx5=shw.vtx.z;
  }

/*
  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,5);
  nu.shwExists6=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists6=true;
    nu.shwEnCor6=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor6=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnMip6=shw.ph.mip;
    nu.planeShwBeg6=shw.plane.beg;
    nu.planeShwEnd6=shw.plane.end;
    nu.xShwVtx6=shw.vtx.x;
    nu.yShwVtx6=shw.vtx.y;
    nu.zShwVtx6=shw.vtx.z;
  }

  pshw=lib.reco.GetShowerWithIndexX(ntp,evt,6);
  nu.shwExists7=false;
  //check if shower exists
  if (pshw) {
    //get a reference
    const NtpSRShower& shw=*pshw;
    //fill the variable(s)
    nu.shwExists7=true;
    nu.shwEnCor7=lib.reco.GetShowerEnergyCor(shw.shwph.linCCgev,CandShowerHandle::kCC,nu);
    nu.shwEnNoCor7=shw.shwph.linCCgev;//shw.ph.gev=shw.shwph.wtCCgev
    nu.shwEnMip7=shw.ph.mip;
    nu.planeShwBeg7=shw.plane.beg;
    nu.planeShwEnd7=shw.plane.end;
    nu.xShwVtx7=shw.vtx.x;
    nu.yShwVtx7=shw.vtx.y;
    nu.zShwVtx7=shw.vtx.z;
  }
*/
}

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

void NuExtraction::ExtractPrimaryTrkInfo(const NtpStRecord& ntp,
                                         const NtpSREvent& evt,
                                         NuEvent& nu) const
{
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //get pointer to primary trk
  const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,0);
  
  //set this to default to no track exists
  nu.trkExists1=false;

  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    nu.trkExists1=true;
    nu.trkIndex1=trk.index;
    nu.ndigitTrk1=trk.ndigit;
    nu.nstripTrk1=trk.nstrip;
    nu.trkEnCorRange1=lib.reco.GetTrackEnergyFromRangeCor
      (trk.momentum.range,nu);
    nu.trkEnCorCurv1=lib.reco.GetTrackEnergyFromCurvCor
      (trk.momentum.qp,nu);
    nu.trkShwEnNear1 = lib.reco.GetShowerEnergyNearTrack(ntp, trk,
                                                         &nu.trkShwEnNearDW1);
    nu.trkMomentumRange1=trk.momentum.range;
    nu.containedTrk1=trk.contained;
    nu.trkfitpass1=trk.fit.pass;
    nu.trkvtxdcosz1=trk.vtx.dcosz;
    nu.trkvtxdcosy1=trk.vtx.dcosy;
    nu.trknplane1=trk.plane.n;
    nu.charge1=lib.reco.GetTrackCharge(trk, nu.hornIsReverse);
    nu.qp1=trk.momentum.qp;
    nu.qp_rangebiased1=trk.momentum.qp_rangebiased;
    nu.sigqp1=trk.momentum.eqp;
    if (nu.sigqp1 != 0.0) nu.qp_sigqp1= nu.qp1 / nu.sigqp1;
    nu.chi21=trk.fit.chi2;
    nu.ndof1=trk.fit.ndof;
    nu.qpFraction1=lib.reco.GetTrkQPFraction(trk);
    nu.trkVtxUVDiffPl1=trk.plane.begu-trk.plane.begv;
    nu.trkLength1=abs(trk.plane.end-trk.plane.beg+1);
    nu.planeTrkNu1=trk.plane.nu;
    nu.planeTrkNv1=trk.plane.nv;
    nu.ntrklike1=trk.plane.ntrklike;
    nu.trkphsigcor1=trk.ph.sigcor;
    nu.trkphsigmap1=trk.ph.sigmap;
    nu.trkIdMC1 = lib.reco.GetTrueTrackId(ntp,trk);
    
    //majority curvature is extracted in its own function

    nu.xTrkVtx1=trk.vtx.x;
    nu.yTrkVtx1=trk.vtx.y;
    nu.zTrkVtx1=trk.vtx.z;
    nu.uTrkVtx1=trk.vtx.u;
    nu.vTrkVtx1=trk.vtx.v;
    nu.planeTrkVtx1=trk.vtx.plane;
    nu.planeTrkBeg1=trk.plane.beg;
    nu.planeTrkBegu1=trk.plane.begu;
    nu.planeTrkBegv1=trk.plane.begv;
    nu.stripTrkBeg1   = lib.reco.GetTrackFirstStrip(ntp, trk);
    nu.stripTrkBegu1  = lib.reco.GetTrackFirstUStrip(ntp, trk);
    nu.stripTrkBegv1  = lib.reco.GetTrackFirstVStrip(ntp, trk);
    nu.stripTrkEnd1   = lib.reco.GetTrackLastStrip(ntp,trk);
    nu.stripTrkEndu1  = lib.reco.GetTrackLastUStrip(ntp,trk);
    nu.stripTrkEndv1  = lib.reco.GetTrackLastVStrip(ntp,trk);
    nu.stripTrkBegIsu1= lib.reco.GetTrackFirstStripIsU(ntp, trk);
    nu.stripTrkEndIsu1= lib.reco.GetTrackLastStripIsU(ntp, trk);
    nu.regionTrkVtx1 = lib.reco.GetRegion(ntp,trk.vtx);
    nu.edgeRegionTrkVtx1 = lib.reco.GetEdgeRegion(ntp,trk.vtx);
    nu.edgeRegionTrkEnd1 = lib.reco.GetEdgeRegion(ntp,trk.end);
    nu.phiTrkVtx1 = lib.reco.GetPhi(trk.vtx);
    nu.phiTrkEnd1 = lib.reco.GetPhi(trk.end);
    
    lib.reco.CalculateParallelStripInset(nu, 1);
    lib.reco.CalculatePerpendicularStripFlag(nu, 1);
    lib.reco.CalculateHorizontalVerticalStripNumber(nu, 1);
    
    nu.xTrkEnd1=trk.end.x;
    nu.yTrkEnd1=trk.end.y;
    nu.zTrkEnd1=trk.end.z;
    nu.uTrkEnd1=trk.end.u;
    nu.vTrkEnd1=trk.end.v;
    nu.planeTrkEnd1=trk.plane.end;
    nu.planeTrkEndu1=trk.plane.endu;
    nu.planeTrkEndv1=trk.plane.endv;
    
    nu.drTrkFidall1=trk.fidall.dr;
    nu.dzTrkFidall1=trk.fidall.dz;
    nu.drTrkFidvtx1=trk.fidvtx.dr;
    nu.dzTrkFidvtx1=trk.fidvtx.dz;
    nu.drTrkFidend1=trk.fidend.dr;
    nu.dzTrkFidend1=trk.fidend.dz;
    nu.traceTrkFidall1 = trk.fidall.trace;
    nu.traceTrkFidvtx1 = trk.fidvtx.trace;
    nu.traceTrkFidend1 = trk.fidend.trace;
    
    nu.cosPrTrkVtx1 = lib.reco.GetCosBetweenPr_Theta(trk);
  }
}

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

void NuExtraction::ExtractSecondTrkInfo(const NtpStRecord& ntp,
                                        const NtpSREvent& evt,
                                        NuEvent& nu) const
{
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //get pointer to trk
  const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,1);
  
  //set this to default to no track exists
  nu.trkExists2=false;

  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    nu.trkExists2=true;
    nu.trkIndex2=trk.index;
    nu.ndigitTrk2=trk.ndigit;
    nu.nstripTrk2=trk.nstrip;
    nu.trkEnCorRange2=lib.reco.GetTrackEnergyFromRangeCor
      (trk.momentum.range,nu);
    nu.trkEnCorCurv2=lib.reco.GetTrackEnergyFromCurvCor
      (trk.momentum.qp,nu);
    nu.trkShwEnNear2 = lib.reco.GetShowerEnergyNearTrack(ntp, trk,
                                                         &nu.trkShwEnNearDW2);
    nu.trkMomentumRange2=trk.momentum.range;
    nu.containedTrk2=trk.contained;
    nu.trkfitpass2=trk.fit.pass;
    nu.trkvtxdcosz2=trk.vtx.dcosz;
    nu.trkvtxdcosy2=trk.vtx.dcosy;
    nu.trknplane2=trk.plane.n;
    nu.charge2=lib.reco.GetTrackCharge(trk, nu.hornIsReverse);
    nu.qp2=trk.momentum.qp;
    nu.qp_rangebiased2=trk.momentum.qp_rangebiased;
    nu.sigqp2=trk.momentum.eqp;
    if (nu.sigqp2 != 0.0) nu.qp_sigqp2=nu.qp2 / nu.sigqp2;
    nu.chi22=trk.fit.chi2;
    nu.ndof2=trk.fit.ndof;
    nu.qpFraction2=lib.reco.GetTrkQPFraction(trk);
    nu.trkVtxUVDiffPl2=trk.plane.begu-trk.plane.begv;
    nu.trkLength2=abs(trk.plane.end-trk.plane.beg+1);
    nu.planeTrkNu2=trk.plane.nu;
    nu.planeTrkNv2=trk.plane.nv;
    nu.ntrklike2=trk.plane.ntrklike;
    nu.trkphsigcor2=trk.ph.sigcor;
    nu.trkphsigmap2=trk.ph.sigmap;
    nu.trkIdMC2 = lib.reco.GetTrueTrackId(ntp,trk);

    //majority curvature is extracted in its own function

    nu.xTrkVtx2=trk.vtx.x;
    nu.yTrkVtx2=trk.vtx.y;
    nu.zTrkVtx2=trk.vtx.z;
    nu.uTrkVtx2=trk.vtx.u;
    nu.vTrkVtx2=trk.vtx.v;
    nu.planeTrkVtx2=trk.vtx.plane;
    nu.planeTrkBeg2=trk.plane.beg;
    nu.planeTrkBegu2=trk.plane.begu;
    nu.planeTrkBegv2=trk.plane.begv;
    nu.stripTrkBeg2   = lib.reco.GetTrackFirstStrip(ntp, trk);
    nu.stripTrkBegu2  = lib.reco.GetTrackFirstUStrip(ntp, trk);
    nu.stripTrkBegv2  = lib.reco.GetTrackFirstVStrip(ntp, trk);
    nu.stripTrkEnd2   = lib.reco.GetTrackLastStrip(ntp,trk);
    nu.stripTrkEndu2  = lib.reco.GetTrackLastUStrip(ntp,trk);
    nu.stripTrkEndv2  = lib.reco.GetTrackLastVStrip(ntp,trk);
    nu.stripTrkBegIsu2= lib.reco.GetTrackFirstStripIsU(ntp, trk);
    nu.stripTrkEndIsu2= lib.reco.GetTrackLastStripIsU(ntp, trk);
    nu.regionTrkVtx2 = lib.reco.GetRegion(ntp,trk.vtx);
    nu.edgeRegionTrkVtx2 = lib.reco.GetEdgeRegion(ntp,trk.vtx);
    nu.edgeRegionTrkEnd2 = lib.reco.GetEdgeRegion(ntp,trk.end);
    nu.phiTrkVtx2 = lib.reco.GetPhi(trk.vtx); //TMath::ATan2(trk.vtx.y, trk.vtx.x);
    nu.phiTrkEnd2 = lib.reco.GetPhi(trk.end);

    lib.reco.CalculateParallelStripInset(nu, 2);
    lib.reco.CalculatePerpendicularStripFlag(nu, 2);
    lib.reco.CalculateHorizontalVerticalStripNumber(nu, 2);
    
    nu.xTrkEnd2=trk.end.x;
    nu.yTrkEnd2=trk.end.y;
    nu.zTrkEnd2=trk.end.z;
    nu.uTrkEnd2=trk.end.u;
    nu.vTrkEnd2=trk.end.v;
    nu.planeTrkEnd2=trk.plane.end;
    nu.planeTrkEndu2=trk.plane.endu;
    nu.planeTrkEndv2=trk.plane.endv;

    nu.drTrkFidall2=trk.fidall.dr;
    nu.dzTrkFidall2=trk.fidall.dz;
    nu.drTrkFidvtx2=trk.fidvtx.dr;
    nu.dzTrkFidvtx2=trk.fidvtx.dz;
    nu.drTrkFidend2=trk.fidend.dr;
    nu.dzTrkFidend2=trk.fidend.dz;
    nu.traceTrkFidall2 = trk.fidall.trace;
    nu.traceTrkFidvtx2 = trk.fidvtx.trace;
    nu.traceTrkFidend2 = trk.fidend.trace;
    
    nu.cosPrTrkVtx2 = lib.reco.GetCosBetweenPr_Theta(trk);
    
  }
}

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

void NuExtraction::ExtractThirdTrkInfo(const NtpStRecord& ntp,
                                       const NtpSREvent& evt,
                                       NuEvent& nu) const
{
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //get pointer to trk
  const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,2);
  
  //set this to default to no track exists
  nu.trkExists3=false;

  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    nu.trkExists3=true;
    nu.trkIndex3=trk.index;
    nu.ndigitTrk3=trk.ndigit;
    nu.nstripTrk3=trk.nstrip;
    nu.trkEnCorRange3=lib.reco.GetTrackEnergyFromRangeCor
      (trk.momentum.range,nu);
    nu.trkEnCorCurv3=lib.reco.GetTrackEnergyFromCurvCor
      (trk.momentum.qp,nu);
    nu.trkShwEnNear3 = lib.reco.GetShowerEnergyNearTrack(ntp, trk,
                                                         &nu.trkShwEnNearDW3);
    nu.trkMomentumRange3=trk.momentum.range;
    nu.containedTrk3=trk.contained;
    nu.trkfitpass3=trk.fit.pass;
    nu.trkvtxdcosz3=trk.vtx.dcosz;
    nu.trkvtxdcosy3=trk.vtx.dcosy;
    nu.trknplane3=trk.plane.n;
    nu.charge3=lib.reco.GetTrackCharge(trk, nu.hornIsReverse);
    nu.qp3=trk.momentum.qp;
    nu.qp_rangebiased3=trk.momentum.qp_rangebiased;
    nu.sigqp3=trk.momentum.eqp;
    if (nu.sigqp3 != 0.0) nu.qp_sigqp3=nu.qp3 / nu.sigqp3;
    nu.chi23=trk.fit.chi2;
    nu.ndof3=trk.fit.ndof;
    nu.qpFraction3=lib.reco.GetTrkQPFraction(trk);
    nu.trkVtxUVDiffPl3=trk.plane.begu-trk.plane.begv;
    nu.trkLength3=abs(trk.plane.end-trk.plane.beg+1);
    nu.planeTrkNu3=trk.plane.nu;
    nu.planeTrkNv3=trk.plane.nv;
    nu.ntrklike3=trk.plane.ntrklike;
    nu.trkphsigcor3=trk.ph.sigcor;
    nu.trkphsigmap3=trk.ph.sigmap;
    nu.trkIdMC3 = lib.reco.GetTrueTrackId(ntp,trk);
    
    //majority curvature is extracted in its own function

    nu.xTrkVtx3=trk.vtx.x;
    nu.yTrkVtx3=trk.vtx.y;
    nu.zTrkVtx3=trk.vtx.z;
    nu.uTrkVtx3=trk.vtx.u;
    nu.vTrkVtx3=trk.vtx.v;
    nu.planeTrkVtx3=trk.vtx.plane;
    nu.planeTrkBeg3=trk.plane.beg;
    nu.planeTrkBegu3=trk.plane.begu;
    nu.planeTrkBegv3=trk.plane.begv;
    nu.stripTrkBeg3   = lib.reco.GetTrackFirstStrip(ntp, trk);
    nu.stripTrkBegu3  = lib.reco.GetTrackFirstUStrip(ntp, trk);
    nu.stripTrkBegv3  = lib.reco.GetTrackFirstVStrip(ntp, trk);
    nu.stripTrkEnd3   = lib.reco.GetTrackLastStrip(ntp,trk);
    nu.stripTrkEndu3  = lib.reco.GetTrackLastUStrip(ntp,trk);
    nu.stripTrkEndv3  = lib.reco.GetTrackLastVStrip(ntp,trk);
    nu.stripTrkBegIsu3= lib.reco.GetTrackFirstStripIsU(ntp, trk);
    nu.stripTrkEndIsu3= lib.reco.GetTrackLastStripIsU(ntp, trk);
    nu.regionTrkVtx3 = lib.reco.GetRegion(ntp, trk.vtx);
    nu.edgeRegionTrkVtx3 = lib.reco.GetEdgeRegion(ntp,trk.vtx);
    nu.edgeRegionTrkEnd3 = lib.reco.GetEdgeRegion(ntp,trk.end);
    nu.phiTrkVtx3 = lib.reco.GetPhi(trk.vtx);
    nu.phiTrkEnd3 = lib.reco.GetPhi(trk.end);

    lib.reco.CalculateParallelStripInset(nu, 3);
    lib.reco.CalculatePerpendicularStripFlag(nu, 3);
    lib.reco.CalculateHorizontalVerticalStripNumber(nu, 3);
    
    nu.xTrkEnd3=trk.end.x;
    nu.yTrkEnd3=trk.end.y;
    nu.zTrkEnd3=trk.end.z;
    nu.uTrkEnd3=trk.end.u;
    nu.vTrkEnd3=trk.end.v;
    nu.planeTrkEnd3=trk.plane.end;
    nu.planeTrkEndu3=trk.plane.endu;
    nu.planeTrkEndv3=trk.plane.endv;

    nu.drTrkFidall3=trk.fidall.dr;
    nu.dzTrkFidall3=trk.fidall.dz;
    nu.drTrkFidvtx3=trk.fidvtx.dr;
    nu.dzTrkFidvtx3=trk.fidvtx.dz;
    nu.drTrkFidend3=trk.fidend.dr;
    nu.dzTrkFidend3=trk.fidend.dz;
    nu.traceTrkFidall3 = trk.fidall.trace;
    nu.traceTrkFidvtx3 = trk.fidvtx.trace;
    nu.traceTrkFidend3 = trk.fidend.trace;
    
    nu.cosPrTrkVtx3 = lib.reco.GetCosBetweenPr_Theta(trk);
  }
}

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

void NuExtraction::ExtractMajorityCurvature(const NtpStRecord& ntp,
                                            const NtpSREvent& evt,
                                            NuEvent& nu) const
{
  //check if a track exists
  if (nu.ntrk<=0) {
    MAXMSG("NuReco",Msg::kInfo,3)
      <<"ExtractMajorityCurvature: no tracks so nothing to do"
      <<", ntrk="<<nu.ntrk<<endl;
    nu.trkExists=false;
    return;
  }

  //check if flag is set to not calc maj curv
  if (!nu.calcMajCurv) {
    MAXMSG("NuReco",Msg::kWarning,1)
      <<"ExtractMajorityCurvature: flag set to calcMajCurv="
      <<nu.calcMajCurv<<endl;
    return;
  }
  
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();
  
  //get pointer to 1st trk
  const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,0);
  
  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    //get an instance of the code library
    const NuLibrary& lib=NuLibrary::Instance();
    const MajCInfo majCInfo=lib.majC.CurvatureImproved(&trk);
    nu.jitter1=majCInfo.jitter;
    nu.jPID1=majCInfo.jPID;
    nu.majC1=majCInfo.majC;
    //nu.majCRatio1=majCInfo.majCRatio;
    //nu.rms1=majCInfo.rms;
    //nu.simpleMajC1=majCInfo.simpleMajC;
    nu.smoothMajC1=majCInfo.smoothMajC;
    //nu.sqJitter1=majCInfo.sqJitter;
    //nu.totWidth1=majCInfo.totWidth;
  }

  //get pointer to 2nd trk
  ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,1);
  
  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    //get an instance of the code library
    const NuLibrary& lib=NuLibrary::Instance();
    const MajCInfo majCInfo=lib.majC.CurvatureImproved(&trk);
    nu.jitter2=majCInfo.jitter;
    nu.jPID2=majCInfo.jPID;
    nu.majC2=majCInfo.majC;
    //nu.majCRatio2=majCInfo.majCRatio;
    //nu.rms2=majCInfo.rms;
    //nu.simpleMajC2=majCInfo.simpleMajC;
    nu.smoothMajC2=majCInfo.smoothMajC;
    //nu.sqJitter2=majCInfo.sqJitter;
    //nu.totWidth2=majCInfo.totWidth;
  }

  //get pointer to 3rd trk
  ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,2);
  
  //check if track exists
  if (ptrk) {
    //get a reference
    const NtpSRTrack& trk=*ptrk;
    
    //get an instance of the code library
    const NuLibrary& lib=NuLibrary::Instance();
    const MajCInfo majCInfo=lib.majC.CurvatureImproved(&trk);
    nu.jitter3=majCInfo.jitter;
    nu.jPID3=majCInfo.jPID;
    nu.majC3=majCInfo.majC;
    //nu.majCRatio3=majCInfo.majCRatio;
    //nu.rms3=majCInfo.rms;
    //nu.simpleMajC3=majCInfo.simpleMajC;
    nu.smoothMajC3=majCInfo.smoothMajC;
    //nu.sqJitter3=majCInfo.sqJitter;
    //nu.totWidth3=majCInfo.totWidth;
  }

  //copy across the MajCurv variables for the best track
  lib.reco.SetBestTrkMajorityCurvature(nu);
}

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

void NuExtraction::ExtractAdditionalDSTVariables(NuEvent& /*nu*/) const
{

/*
  //get an instance of the code library
  const NuLibrary& lib=NuLibrary::Instance();

  //set fiducial volume variables useful for DST interactive use  
  if (nu.detector==Detector::kNear) {
    nu.fvnmb=lib.cuts.IsInFidVolNDNuMuBar
      (nu.xTrkVtx,nu.yTrkVtx,nu.zTrkVtx);
    nu.fvpitt=lib.cuts.IsInFidVolPitt
      (nu.xTrkVtx,nu.yTrkVtx,nu.zTrkVtx,nu.uTrkVtx,nu.vTrkVtx);
    nu.fvcc=lib.cuts.IsInFidVolNDCC0250Std
      (nu.xTrkVtx,nu.yTrkVtx,nu.zTrkVtx);
  }
  else if (nu.detector==Detector::kFar) {
    nu.fvcc=lib.cuts.IsInFidVolFDCC0250Std
      (nu.planeTrkVtx,nu.rTrkVtx);
    //just set these other two to the std CC one
    //there are no specialised FD fid vols yet
    nu.fvnmb=nu.fvcc;
    nu.fvpitt=nu.fvcc;
  }
  else {
    cout<<"Ahhhh, det="<<nu.detector<<endl;
  }
*/

}

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

void NuExtraction::ExtractMinMaxEvtTimes(const NtpStRecord& ntp,
                                         const NtpSREvent& evt,
                                         NuEvent& nu) const
{
  //get a multiset to hold all the strip times
  multiset<Double_t> times;
  const TClonesArray& stpTca=(*ntp.stp);

  Double_t timemax=-1.e10;
  Double_t timemin=1.e10;

  MsgFormat ffmt("%9.9f");

  //Extract the crate t0:

  NtpSRTimeStatus timeStatus = ntp.timestatus;
  nu.crateT0 = timeStatus.crate_t0_ns * Munits::ns;
  //  cout << "Extracted crate t0: " << nu.crateT0 << endl;

/*
  Double_t snarlTime=(nu.timeNanoSec*(Munits::ns));
  if (snarlTime>0.999999981 && snarlTime<0.999999983) {
    const Int_t numStps=stpTca.GetEntriesFast();
    Int_t counter=0;
    //loop over strips in snarl
    for (Int_t i=0;i<numStps;++i) {
      const NtpSRStrip* pstp=
        dynamic_cast<NtpSRStrip*>(stpTca[i]);
      const NtpSRStrip& stp=(*pstp);
      
      if (stp.time1>0.000029999) {
        MAXMSG("NuExtraction",Msg::kInfo,10000)
          <<"run="<<nu.run<<", subrun="<<nu.subRun
          <<", snarl="<<nu.snarl<<", i="<<i
          <<", count="<<counter<<", evt="<<nu.evt
          <<", ind="<<stp.index<<", s,pl="<<stp.strip<<","<<stp.plane
          <<endl
          <<"  t1="<<ffmt(stp.time1)
          <<", snlTime="<<ffmt(snarlTime)
          <<", raw="<<stp.ph1.raw
          <<", sigC="<<stp.ph1.sigcor
          <<endl;
      }
      counter++;
    }
  }
*/

  MAXMSG("NuExtraction",Msg::kDebug,200)
    <<"evt.nstrip="<<evt.nstrip<<endl;
  for (Int_t i=0;i<evt.nstrip;i++) {        
    //check for bug where strip index is -1
    if (evt.stp[i]<0) {
      MAXMSG("MeuCuts",Msg::kInfo,50)
        <<"Skipping strip with evt.stp[i]="<<evt.stp[i]<<endl;
      continue;
    }
    
    const NtpSRStrip& stp=
      *(dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]));

    Double_t snarlTime=(nu.timeNanoSec*(Munits::ns));
    Double_t striptime=-999;

    //get the min/max time
    //code from MadTVAnalysis
    if (nu.detector==Detector::kNear) {
      if (nu.simFlag==SimFlag::kMC &&
          snarlTime>0.999999981 && snarlTime<0.999999983) {
        MAXMSG("NuExtraction",Msg::kWarning,1)
          <<"Deal with odd times: setting fNanoSec to zero, snarlTime="
          <<ffmt(snarlTime)<<endl;
        snarlTime=0;
      }
      
      //get the time relative to the snarl time
      striptime=stp.time1-snarlTime;

      //check if min/max
      if(striptime<=timemin && striptime!=-999) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
    else if (nu.detector==Detector::kFar) {
      //get the times relative to the snarl time
      Double_t striptime1=stp.time1-snarlTime;
      Double_t striptime0=stp.time0-snarlTime;
      striptime=-999;
      if(striptime1>0 && striptime0<0) striptime=striptime1;
      if(striptime0>0 && striptime1<0) striptime=striptime0;
      if(striptime0>0 && striptime1>0) striptime=(striptime0+
                                                  striptime1)/2.;

      //check if min/max
      if(striptime<=timemin && striptime!=-999) timemin=striptime;
      if(striptime>=timemax) timemax=striptime;
    }
    else cout<<"Ahhh, detector not known"<<endl;

    Double_t time=stp.time0;
    if (time<0 || time>1) time=stp.time1;
    
    if (time>0 && time<=1) {
      times.insert(time);
    }
    //else just don't put the time in the map - clearly rubbish
    if (nu.detector==Detector::kNear) {
      MAXMSG("NuExtraction",Msg::kDebug,1000)
        <<"i="<<stp.index<<", s="<<stp.strip<<", t1="<<ffmt(stp.time1)
        <<", tmn="<<ffmt(timemin)<<", tmx="<<ffmt(timemax)
        <<", snlTime="<<ffmt(snarlTime)<<endl;
    }
    else {
      MAXMSG("NuExtraction",Msg::kDebug,1000)
        <<"i="<<stp.index<<", s="<<stp.strip
        <<", t0="<<ffmt(stp.time0)<<", t1="<<ffmt(stp.time1)
        <<", tmn="<<ffmt(timemin)<<", tmx="<<ffmt(timemax)
        <<", snlTime="<<ffmt(snarlTime)
        //<<", pl="<<stp.plane
        //<<", sigC="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
        <<endl;
    }
  }
  
  //get the median time from the map
  multiset<Double_t>::const_iterator it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTime=*it;
  MAXMSG("NuExtraction",Msg::kDebug,100)
    <<"Median time="<<medianTime<<endl;
  
  //set the NuEvent variables
  nu.medianTime=medianTime;
  nu.timeEvtMin=timemin;
  nu.timeEvtMax=timemax;


  //sanity check and debug printout
  if ((timemin+nu.crateT0)<0 || (nu.run==36608 && nu.snarl==42717 && false)) {
    //Bad time, run=36608, snarl=42717, evt=0
    Double_t timemax=-1.e-10;
    Double_t timemin=1.e10;

    MAXMSG("NuExtraction",Msg::kError,1000)
      <<"Bad time, run="<<nu.run<<", snarl="<<nu.snarl
      <<", evt="<<nu.evt<<endl;
    for (Int_t i=0;i<evt.nstrip;i++) {
      //check for bug where strip index is -1
      if (evt.stp[i]<0) {
        MAXMSG("MeuCuts",Msg::kInfo,50)
          <<"Skipping strip with evt.stp[i]="<<evt.stp[i]<<endl;
        continue;
      }
      
      const NtpSRStrip& stp=
        *(dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]));
      
      Double_t snarlTime=(nu.timeNanoSec*(Munits::ns));
      Double_t striptime=-999;
      
      //get the min/max time
      //code from MadTVAnalysis
      if (nu.detector==Detector::kNear) {
        //get the time relative to the snarl time
        striptime=stp.time1-snarlTime;
        
        //check if min/max
        if(striptime<=timemin) timemin=striptime;
        if(striptime>=timemax) timemax=striptime;
      }
      else if (nu.detector==Detector::kFar) {
        //get the times relative to the snarl time
        Double_t striptime1=stp.time1-snarlTime;
        Double_t striptime0=stp.time0-snarlTime;
        striptime=-999;
        if(striptime1>0 && striptime0<0) striptime=striptime1;
        if(striptime0>0 && striptime1<0) striptime=striptime0;
        if(striptime0>0 && striptime1>0) striptime=(striptime0+
                                                    striptime1)/2.;
        
        //check if min/max
        if(striptime<=timemin) timemin=striptime;
        if(striptime>=timemax) timemax=striptime;
      }
      else cout<<"Ahhh, detector not known"<<endl;
      
      MAXMSG("NuExtraction",Msg::kInfo,1000)
        <<"strip="<<stp.strip
        <<", t0="<<ffmt(stp.time0)<<", t1="<<ffmt(stp.time1)
        <<", stp.index="<<stp.index
        <<endl
        <<"  tmn="<<ffmt(timemin)<<", tmx="<<ffmt(timemax)
        <<", stime="<<ffmt(striptime)
        <<", snlTime="<<ffmt(snarlTime)<<endl
        <<"  plane="<<stp.plane
        <<", ph0.raw="<<stp.ph0.raw
        <<", ph1.raw="<<stp.ph1.raw
        <<endl;
    }
  }
}

//......................................................................
void NuExtraction::ExtractNCInfo(const NtpStRecord* pntp,
                   Int_t eventNb,
                   NuEvent& nu) const
{
  // Protection against reset of gDirectory
  // The call to InitializekNN() uses another file.
  TDirectory *tmpd = gDirectory;
  MSG("NuExtraction",Msg::kDebug) << "Start ExtractNCInfo, gDirectory is : ";
  //gDirectory->pwd();

  // See AnalysisNtuple/Module/CondensedNtpModuleNC::Ana()
  // create the ntpManipulator object 
  static ANtpInfoObjectFillerNC   *InfoFiller   = new ANtpInfoObjectFillerNC();
  static ANtpRecoNtpManipulator *ntpManipulator = new ANtpRecoNtpManipulator();
  NtpStRecord            *tryRec = const_cast<NtpStRecord*>(pntp);
  ntpManipulator->SetRecord(tryRec);
  static ANtpEventInfoNC   *EventInfo  = new ANtpEventInfoNC();
  static ANtpTrackInfoNC   *TrackInfo  = new ANtpTrackInfoNC();
  static ANtpShowerInfoNC  *ShowerInfo = new ANtpShowerInfoNC();
  static ANtpTruthInfoBeam *TruthInfo  = new ANtpTruthInfoBeam();                                                                                             

  // Extract Header Info
  //====================
  // Detector - 1 near, 2 far - needed for filling the HeaderInfo
  Detector::Detector_t det = tryRec->GetVldContext()->GetDetector();
  InfoFiller->SetDetector(det);
  InfoFiller->SetStripArray(ntpManipulator->GetStripArray());
  InfoFiller->SetClusterArray(ntpManipulator->GetClusterArray());
  
  // Extract Event/Track/Shower/Truth Info
  //======================================
  //set up which flags you want to use to determine the primary shower or track
  //a value of 0 for a flag means it will not be used
  ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height
  ntpManipulator->SetPrimaryShowerCriteria(0,1);  // nplanes, total pulse height
                                                                                             

  // Initialize kNN, only once per snarl
  //====================================
  InfoFiller->InitializekNN(ntpManipulator);

  // Reset the Info Objects
  //=======================
  EventInfo->Reset();
  TrackInfo->Reset();
  ShowerInfo->Reset();
  TruthInfo->Reset();                                                                                             

  // Fill event, track, shower and truth and calculate variables needed
  //===================================================================
  InfoFiller->FillInformation(eventNb, ntpManipulator, EventInfo, TrackInfo, ShowerInfo, TruthInfo);
  
  // Fill the NuEvent variables
  //===========================
  // event
  nu.closeTimeDeltaZ = EventInfo->closeTimeDeltaZ;
  nu.edgeActivityStrips = EventInfo->edgeActivityStrips;
  nu.edgeActivityPH = EventInfo->edgeActivityPH;
  nu.oppEdgeStrips = EventInfo->oppEdgeStrips;
  nu.oppEdgePH = EventInfo->oppEdgePH;
  nu.vtxMetersToCoilEvt = EventInfo->vtxMetersToCoil;
  nu.vtxMetersToCloseEdgeEvt = EventInfo->vtxMetersToCloseEdge;
  nu.minTimeSeparation = EventInfo->minTimeSeparation;                                                                                             

  // shower
  nu.transverseRMSU = ShowerInfo->transverseRMSU;
  nu.transverseRMSV = ShowerInfo->transverseRMSV;

  // track
  nu.dtdz = TrackInfo->dtdz;
  nu.endMetersToCloseEdge = TrackInfo->endMetersToCloseEdge;
  nu.vtxMetersToCloseEdgeTrk = TrackInfo->vtxMetersToCloseEdge;
  nu.vtxMetersToCoilTrk = TrackInfo->vtxMetersToCoil;
  nu.traceEndZ = TrackInfo->traceEndZ;                                                                                             
    
  MSG("NuExtraction",Msg::kDebug) << "End ExtractNCInfo, gDirectory is : ";
  gDirectory = tmpd;
  //gDirectory->pwd();
}

//......................................................................
void NuExtraction::NuMCEventFromNuEvent(const NuEvent& nu,
                                        NuMCEvent& numc) const
{
  //copy the info across

  //book keeping  quantities
  numc.entry=nu.entry;

  //snarl/run based quantities
  numc.run=nu.run;
  numc.subRun=nu.subRun;
  numc.snarl=nu.snarl;    

  numc.runPeriod=nu.runPeriod;
  numc.intensity=nu.intensity;
  numc.hornIsReverse=nu.hornIsReverse;
  numc.beamType=nu.beamType;

  numc.detector=nu.detector;
  numc.simFlag=nu.simFlag;
  numc.timeSec=nu.timeSec;
  numc.timeNanoSec=nu.timeNanoSec;
  numc.timeSeconds=nu.timeSeconds;
  numc.trigtime=nu.trigtime;

  numc.releaseType=nu.releaseType;
  numc.anaVersion=nu.anaVersion;

  //truth variables
  numc.energyMC=nu.energyMC;

  numc.neuEnMC=nu.neuEnMC;
  numc.neuPxMC=nu.neuPxMC;
  numc.neuPyMC=nu.neuPyMC;
  numc.neuPzMC=nu.neuPzMC;

  numc.mu1EnMC=nu.mu1EnMC;
  numc.mu1PxMC=nu.mu1PxMC;
  numc.mu1PyMC=nu.mu1PyMC;
  numc.mu1PzMC=nu.mu1PzMC;

  numc.tgtEnMC=nu.tgtEnMC;
  numc.tgtPxMC=nu.tgtPxMC;
  numc.tgtPyMC=nu.tgtPyMC;
  numc.tgtPzMC=nu.tgtPzMC;

  numc.zMC=nu.zMC;
  numc.aMC=nu.aMC;
  numc.nucleusMC=nu.nucleusMC;
  numc.initialStateMC=nu.initialStateMC;
  numc.hadronicFinalStateMC=nu.hadronicFinalStateMC;
  numc.numPreInukeFSprotMC=nu.numPreInukeFSprotMC;
  numc.numPreInukeFSneutMC=nu.numPreInukeFSneutMC;
  numc.maxMomPreInukeFSprotMC=numc.maxMomPreInukeFSprotMC;
  numc.maxMomPreInukeFSneutMC=numc.maxMomPreInukeFSneutMC;

  numc.yMC=nu.yMC;
  numc.y2MC=nu.y2MC;
  numc.xMC=nu.xMC;
  numc.q2MC=nu.q2MC;
  numc.w2MC=nu.w2MC;

  numc.trkEnMC=nu.trkEnMC;
  numc.trkEn2MC=nu.trkEn2MC;
  numc.shwEnMC=nu.shwEnMC;
  numc.shwEn2MC=nu.shwEn2MC;

  numc.trkEndEnMC=nu.trkEndEnMC;
  numc.trkStartEnMC=nu.trkStartEnMC;
  numc.trkContainmentMC=nu.trkContainmentMC;

  numc.sigma=nu.sigma;
  numc.iaction=nu.iaction;
  numc.iresonance=nu.iresonance;
  numc.inu=nu.inu;
  numc.inunoosc=nu.inunoosc;
  numc.itg=nu.itg;
    
  numc.vtxxMC=nu.vtxxMC;
  numc.vtxyMC=nu.vtxyMC;
  numc.vtxzMC=nu.vtxzMC;
  numc.vtxuMC=nu.vtxuMC;
  numc.vtxvMC=nu.vtxvMC;
  numc.planeTrkVtxMC=nu.planeTrkVtxMC;
  numc.rTrkVtxMC=nu.rTrkVtxMC;

  numc.mc=nu.mc;//copy across the index of the true event

  numc.Npz=nu.Npz;
  numc.NdxdzNea=nu.NdxdzNea;
  numc.NdydzNea=nu.NdydzNea;
  numc.NenergyN=nu.NenergyN;
  numc.NWtNear=nu.NWtNear;
  numc.NdxdzFar=nu.NdxdzFar;
  numc.NdydzFar=nu.NdydzFar;
  numc.NenergyF=nu.NenergyF;
  numc.NWtFar=nu.NWtFar;
  numc.Ndecay=nu.Ndecay;
  numc.Vx=nu.Vx;
  numc.Vy=nu.Vy;
  numc.Vz=nu.Vz;
  numc.pdPx=nu.pdPx;
  numc.pdPy=nu.pdPy;
  numc.pdPz=nu.pdPz;
  numc.ppdxdz=nu.ppdxdz;
  numc.ppdydz=nu.ppdydz;
  numc.pppz=nu.pppz;
  numc.ppenergy=nu.ppenergy;
  numc.ppmedium=nu.ppmedium;
  numc.ppvx=nu.ppvx;
  numc.ppvy=nu.ppvy;
  numc.ppvz=nu.ppvz;
  numc.ptype=nu.ptype;
  numc.Necm=nu.Necm;
  numc.Nimpwt=nu.Nimpwt;
  numc.tvx=nu.tvx;
  numc.tvy=nu.tvy;
  numc.tvz=nu.tvz;
  numc.tpx=nu.tpx;
  numc.tpy=nu.tpy;
  numc.tpz=nu.tpz;
  numc.tptype=nu.tptype;
  numc.tgen=nu.tgen;

  //copy the weights
  numc.rw=nu.rw;
  numc.fluxErr=nu.fluxErr;
  numc.rwActual=nu.rwActual;
  numc.generatorWeight=nu.generatorWeight;
  numc.detectorWeight=nu.detectorWeight;

  numc.trkEnWeight=nu.trkEnWeight;
  numc.shwEnWeight=nu.shwEnWeight;
  numc.beamWeight=nu.beamWeight;
  numc.fluxErrHadProdAfterTune=nu.fluxErrHadProdAfterTune;
  numc.fluxErrTotalErrorPreTune=nu.fluxErrTotalErrorPreTune;
  numc.fluxErrTotalErrorAfterTune=nu.fluxErrTotalErrorAfterTune;
  numc.detectorWeightNMB=nu.detectorWeightNMB;
  numc.detectorWeightNM=nu.detectorWeightNM;

  numc.trkEnWeightRunI=nu.trkEnWeightRunI;
  numc.shwEnWeightRunI=nu.shwEnWeightRunI;
  numc.beamWeightRunI=nu.beamWeightRunI;
  numc.fluxErrHadProdAfterTuneRunI=nu.fluxErrHadProdAfterTuneRunI;
  numc.fluxErrTotalErrorPreTuneRunI=nu.fluxErrTotalErrorPreTuneRunI;
  numc.fluxErrTotalErrorAfterTuneRunI=nu.fluxErrTotalErrorAfterTuneRunI;
  numc.detectorWeightNMBRunI=nu.detectorWeightNMBRunI;
  numc.detectorWeightNMRunI=nu.detectorWeightNMRunI;

  numc.trkEnWeightRunII=nu.trkEnWeightRunII;
  numc.shwEnWeightRunII=nu.shwEnWeightRunII;
  numc.beamWeightRunII=nu.beamWeightRunII;
  numc.fluxErrHadProdAfterTuneRunII=nu.fluxErrHadProdAfterTuneRunII;
  numc.fluxErrTotalErrorPreTuneRunII=nu.fluxErrTotalErrorPreTuneRunII;
  numc.fluxErrTotalErrorAfterTuneRunII=
    nu.fluxErrTotalErrorAfterTuneRunII;
  numc.detectorWeightNMBRunII=nu.detectorWeightNMBRunII;
  numc.detectorWeightNMRunII=nu.detectorWeightNMRunII;
  
  numc.InukeNwts  =  nu.InukeNwts ;
  numc.InukePiCExchgP  =  nu.InukePiCExchgP ;   //0  
  numc.InukePiCExchgN  =  nu.InukePiCExchgN ;   //1  
  numc.InukePiEScatP  =  nu.InukePiEScatP ;     //2  
  numc.InukePiEScatN  =  nu.InukePiEScatN ;     //3  
  numc.InukePiInEScatP  =  nu.InukePiInEScatP ;  //4  
  numc.InukePiInEScatN  =  nu.InukePiInEScatN ;  //5  
  numc.InukePiAbsorbP  =  nu.InukePiAbsorbP ;   //6  
  numc.InukePiAbsorbN  =  nu.InukePiAbsorbN ;   //7  
  numc.InukePi2PiP  =  nu.InukePi2PiP ;         //8  
  numc.InukePi2PiN  =  nu.InukePi2PiN ;         //9  
  numc.InukeNknockP  =  nu.InukeNknockP ;       //10 
  numc.InukeNknockN  =  nu.InukeNknockN ;       //11 
  numc.InukeNNPiP  =  nu.InukeNNPiP ;           //12 
  numc.InukeNNPiN  =  nu.InukeNNPiN ;           //13 
  numc.InukeFormTP  =  nu.InukeFormTP ;      //14 
  numc.InukeFormTN  =  nu.InukeFormTN ;      //15 
  numc.InukePiXsecP  =  nu.InukePiXsecP ;     //16 
  numc.InukePiXsecN  =  nu.InukePiXsecN ;     //17 
  numc.InukeNXsecP  =  nu.InukeNXsecP ;      //18 
  numc.InukeNXsecN  =  nu.InukeNXsecN ;      //19 
  numc.InukeNucrad  =  nu.InukeNucrad ;
  numc.InukeWrad  =  nu.InukeWrad ;
}

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

---