ANtpInfoObjectFillerNC.cxx

Go to the documentation of this file.

//$Id: ANtpInfoObjectFillerNC.cxx,v 1.47 2009/11/17 18:30:19 tinti Exp $
//
//ANtpInfoObjectFiller.cxx
//
//B. Rebel 12/2005

#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "AnalysisNtuples/Module/ANtpInfoObjectFillerNC.h"
#include "AnalysisNtuples/ANtpTrackInfoNC.h"
#include "AnalysisNtuples/ANtpShowerInfoNC.h"
#include "AnalysisNtuples/ANtpEventInfoNC.h"
#include "MessageService/MsgService.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpDefaultValue.h"
#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRCluster.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "BField/BfieldCoilCurrent.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"
#include "SpillTiming/SpillTimeFinder.h"
#include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"

// #include "AnalysisNtuples/Module/VHS.h"
// #include "AnalysisNtuples/Module/VHSevent.h"

#include "Conventions/Munits.h"
#include "Conventions/PlaneView.h"
#include "Conventions/ReleaseType.h"

#include "Registry/Registry.h"

//ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TSystem.h"

//C++/STL includes
#include <cassert>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <map>
#include <set>

static int kNearPartial = 1;
static int kNearFull = 2;
static int kComplete = 3;
static int kU = 2;
static int kV = 3;

ClassImp(ANtpInfoObjectFillerNC)

CVSID("$Id: ANtpInfoObjectFillerNC.cxx,v 1.47 2009/11/17 18:30:19 tinti Exp $");

//-------------------------------------------------------------------
ANtpInfoObjectFillerNC::ANtpInfoObjectFillerNC() :
  // ANtpInfoObjectFillerBeam(Detector::kUnknown),
  ANtpInfoObjectFillerBeam(Detector::kNear),
  fClusterArray(0),
  fkNNSet(false),
  fVHSPlanes(20),
  fVHSStrips(20)
{
  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "ANtpInfoObjectFillerNC::Default Constructor" << endl;

  FillStripToPixelMaps();

}

//----------------------------------------------------------------------
ANtpInfoObjectFillerNC::ANtpInfoObjectFillerNC(Detector::Detector_t detector) :
  ANtpInfoObjectFillerBeam(detector),
  fClusterArray(0),
  fkNNSet(false),
  fVHSPlanes(20),
  fVHSStrips(20)
{
  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "ANtpInfoObjectFillerNC::Constructor" << endl;

  FillStripToPixelMaps();

}

//----------------------------------------------------------------------
ANtpInfoObjectFillerNC::~ANtpInfoObjectFillerNC()
{
  if(fAnpInterface) delete fAnpInterface;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::FillStripToPixelMaps()
{
  // Fills the strip to pixel map. Call this when you pick a detector

  if(fDetector  == Detector::kFar){

    //Define the FD strip to pixel maps
    Int_t stripToPixelWest[] = {0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6};

    Int_t stripToPixelEast[] = {0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
                                2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0,
                                5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2,
                                7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5,
                                8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7,
                                10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8,
                                13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8, 10,
                                15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8, 10, 13};

    // Fill the strip to pixe map data members
    for(Int_t i = 0; i < 192; ++i){
      fStripToPixelMapWest[i] = stripToPixelWest[i];
      fStripToPixelMapEast[i] = stripToPixelEast[i];
    }

    // set plane coverages
    for(Int_t i = 0; i < 486; ++i) fPlaneCoverage[i] = kComplete;
  } // end filling strip to pixel map for far detector

  else if(fDetector == Detector::kNear){

    // Define the pixel to strip maps
    Int_t pixelToStripU1[] = {40, 53, 22, 35, 11, 41, 54, 23,
                              36, 10, 42, 55, 24, 37,  9, 43,
                              56, 25, 38,  8, 44, 57, 26, 39,
                              7, 45, 58, 27, 19,  6, 46, 59,
                              28, 18,  5, 47, 60, 29, 17,  4,
                              48, 61, 30, 16,  3, 49, 62, 31,
                              15,  2, 50, 63, 32, 14,  1, 51,
                              20, 33, 13,  0, 52, 21, 34, 12};

    Int_t pixelToStripV1[] = {27, 14, 45, 32, 56, 26, 13, 44,
                              31, 57, 25, 12, 43, 30, 58, 24,
                              11, 42, 29, 59, 23, 10, 41, 28,
                              60, 22,  9, 40, 48, 61, 21,  8,
                              39, 49, 62, 20,  7, 38, 50, 63,
                              19,  6, 37, 51, 64, 18,  5, 36,
                              52, 65, 17,  4, 35, 53, 66, 16,
                              47, 34, 54, 67, 15, 46, 33, 55};

    Int_t pixelToStripU2[] = {68, 81, 94, 59, -1, 69, 82, 95,
                              60, -1, 70, 83, 48, 61, -1, 71,
                              84, 49, 62, -1, 72, 85, 50, 63,
                              -1, 73, 86, 51, 64, -1, 74, 87,
                              52, 65, -1, 75, 88, 53, 66, -1,
                              76, 89, 54, 67, -1, 77, 90, 55,
                              -1, -1, 78, 91, 56, -1, -1, 79,
                              92, 57, -1, -1, 80, 93, 58, -1};

    Int_t pixelToStripV2[] = {27, 14,  1, 36, -1, 26, 13,  0,
                              35, -1, 25, 12, 47, 34, -1, 24,
                              11, 46, 33, -1, 23, 10, 45, 32,
                              -1, 22,  9, 44, 31, -1, 21,  8,
                              43, 30, -1, 20,  7, 42, 29, -1,
                              19,  6, 41, 28, -1, 18,  5, 40,
                              -1, -1, 17,  4, 39, -1, -1, 16,
                              3, 38, -1, -1, 15,  2, 37, -1};

    Int_t pixelToStripU3[] = {47, 34, 21,  8, -1, 46, 33, 20,
                              7, -1, 45, 32, 19,  6, -1, 44,
                              31, 18,  5, -1, 43, 30, 17,  4,
                              -1, 42, 29, 16,  3, -1, 41, 28,
                              15,  2, -1, 40, 27, 14,  1, -1,
                              39, 26, 13,  0, -1, 38, 25, 12,
                              -1, -1, 37, 24, 11, -1, -1, 36,
                              23, 10, -1, -1, 35, 22,  9, -1};

    Int_t pixelToStripV3[] = {48, 61, 74, 87, -1, 49, 62, 75,
                              88, -1, 50, 63, 76, 89, -1, 51,
                              64, 77, 90, -1, 52, 65, 78, 91,
                              -1, 53, 66, 79, 92, -1, 54, 67,
                              80, 93, -1, 55, 68, 81, 94, -1,
                              56, 69, 82, 95, -1, 57, 70, 83,
                              -1, -1, 58, 71, 84, -1, -1, 59,
                              72, 85, -1, -1, 60, 73, 86, -1};

    //now fill the strip to pixel map data members
    for(Int_t i = 0; i < 64; ++i){
      fStripToPixelMapNearU1[pixelToStripU1[i]] = i;
      fStripToPixelMapNearV1[pixelToStripV1[i]] = i;

      fStripToPixelMapNearU2[i] = -1;
      if(pixelToStripU2[i]>-1) fStripToPixelMapNearU2[pixelToStripU2[i]] = i;

      fStripToPixelMapNearV2[i] = -1;
      if(pixelToStripV2[i]>-1) fStripToPixelMapNearV2[pixelToStripV2[i]] = i;

      fStripToPixelMapNearU3[i] = -1;
      if(pixelToStripU3[i]>-1) fStripToPixelMapNearU3[pixelToStripU3[i]] = i;

      fStripToPixelMapNearV3[i] = -1;
      if(pixelToStripV3[i]>-1) fStripToPixelMapNearV3[pixelToStripV3[i]] = i;

    }//end fill strip to pixel maps

    //fill the array of plane coverages
    for(Int_t i = 0; i < 282; ++i){
      if( i%5 == 1 ) fPlaneCoverage[i] = kNearFull;
      else fPlaneCoverage[i] = kNearPartial;
    }
  }//end filling strip to pixel map for near detector

  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "ANtpInfoObjectFillerNC::Constructor" << endl;

}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::SetDetector(Detector::Detector_t detector)
{
  // Set the detector. Then fill the strip to pixel maps accordingly
  fDetector = detector;
  FillStripToPixelMaps();
  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::SetClusterArray(const TClonesArray *clusters) {
  fClusterArray = clusters;
}

//----------------------------------------------------------------------
bool ANtpInfoObjectFillerNC::FillInformation(int event,
                                             ANtpRecoNtpManipulator *ntpManipulator,
                                             ANtpEventInfoNC *eventInfo,
                                             ANtpTrackInfoNC *trackInfo,
                                             ANtpShowerInfoNC *showerInfo,
                                             ANtpTruthInfoBeam *truthInfo)
{

  NtpSRTrack  *ntpTrack = 0;
  NtpSREvent  *ntpEvent = 0;
  NtpSRShower *ntpShower = 0;
  NtpMCTruth  *ntpMCTruth = 0;
  NtpTHEvent  *ntpTHEvent = 0;

  MAXMSG("ANtpInfoObjectFillerNC", Msg::kDebug,20)
    << "Begin Fill Information method..." << endl;

  //get event.  the call to ANtpEventManipulator::SetEventInSnarl sets the
  //NtpSREvent data member in that object to the current event so that
  //you can later call for the primary track, shower, etc.
  ntpManipulator->GetEventManipulator()->SetEventInSnarl(event);
  ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent();
  if(!ntpEvent) return false;

  FillEventInformation(ntpManipulator, ntpEvent, eventInfo);
  FillEventTimingAndActivityInformation(ntpManipulator, event, eventInfo);

  //get the primary shower for the event - if no shower is present it
  //returns 0
  ntpShower = ntpManipulator->GetEventManipulator()->GetPrimaryShower();
  if(ntpShower) FillShowerInformation(ntpShower, showerInfo, eventInfo);

  //get the primary track for the event - if no track is present it
  //returns 0
  ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack();
  if(ntpTrack)
    FillTrackInformation(ntpManipulator, ntpTrack, ntpEvent, trackInfo, eventInfo);

  FillCrossOverInformation(ntpTrack, ntpShower, ntpEvent, ntpManipulator,
                           trackInfo, showerInfo, eventInfo);

  // Get the truth information for this track if it is MC
  // Get best neu match from Truthhelper
  ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(event);

  if(ntpTHEvent)
    ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc);
  if(ntpMCTruth)
    FillMCInformation(ntpMCTruth, ntpManipulator->GetStdHepArray(), truthInfo);
  // added by Tobi
  if (truthInfo && ntpTHEvent)
    truthInfo->eventCompleteness = ntpTHEvent->completeall;

  return true;
}


//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::
GetEvtVtxWithFixup(ANtpRecoNtpManipulator* ntpManipulator,
                   int& vtxPlane, float& vtxZ) const
{
  ANtpEventManipulator* evtManip = ntpManipulator->GetEventManipulator();
  const NtpSREvent* ntpEvent = evtManip->GetEvent();

  vtxPlane = ntpEvent->vtx.plane;
  vtxZ = ntpEvent->vtx.z;

  //Correct cedar vertex bug
  const TString reco = ntpManipulator->GetReleaseMCType();
  const bool isCedar = reco.Contains("Cedar");
  if(isCedar){
    NtpVtxFinder::NtpVtxFinder vtxf;
    vtxf.SetTargetEvent(ntpEvent->index, ntpManipulator->GetNtpStRecord());
    if(vtxf.FindVertex() > 0){ // Take fixed vertex if possible
      vtxPlane = vtxf.VtxPlane();
      vtxZ = vtxf.VtxZ();
    }
  } // end if cedar

  const NtpSRShower* ntpShower = evtManip->GetPrimaryShower();
  const NtpSRTrack* ntpTrack = evtManip->GetPrimaryTrack();

  if(ntpTrack && ntpShower){
    if(ntpEvent->ntrack > 0 &&
       ntpTrack->plane.n > ntpShower->plane.n){
      vtxPlane = ntpTrack->vtx.plane;
      vtxZ = ntpTrack->vtx.z - 3.92 * Munits::cm;
    }
  }
}


//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::
GetStripEventTime(const NtpSREvent* ntpEvent,
                  int vtxPlane, float vtxZ,
                  double& evtTime,
                  double& ToFCorrectedEvtTime) const
{
  vector<double> stpTimes;

  ToFCorrectedEvtTime = 0;
  Int_t tfcount=0;
  // loop over strips in current event
  for(int ns = 0; ns < ntpEvent->nstrip; ++ns){
    const int stpIdx = ntpEvent->stp[ns];
    if(stpIdx < 0) continue;
    const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(stpIdx);

    if(ntpStrip->ph1.pe >=2) {
      ToFCorrectedEvtTime += ntpStrip->time1-
        (fabs(ntpStrip->z - vtxZ)/Munits::c_light);
      tfcount++;
    }
    // 5 plane window
    if(ntpStrip->plane >= vtxPlane &&
       ntpStrip->plane < vtxPlane + 5){
      stpTimes.push_back(ntpStrip->time1);
    }
  } // end for ns

   if(tfcount) ToFCorrectedEvtTime /= (double)tfcount;
   else ToFCorrectedEvtTime=ANtpDefVal::kDouble;

  // calculate strip time median for this event
  sort(stpTimes.begin(), stpTimes.end());

  if(stpTimes.size() % 2)
    evtTime = *(stpTimes.begin()+stpTimes.size()/2);
  else
    evtTime = ( *(stpTimes.begin()+stpTimes.size()/2)
                + *(stpTimes.begin()+stpTimes.size()/2-1) )/2.;
}


//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::
FillEventTimingAndActivityInformation(ANtpRecoNtpManipulator* ntpManipulator,
                                      int event,
                                      ANtpEventInfoNC* eventInfo)
{
  ANtpEventManipulator* evtManip = ntpManipulator->GetEventManipulator();

  //the event number passed into the method is the current event we care about
  evtManip->SetEventInSnarl(event);

  const NtpSREvent* ntpEvent = evtManip->GetEvent();

  // VERTEX TIMING AND DELTA Z CALCULATIONS

  int thisVtxPlane; float thisVtxZ;
  GetEvtVtxWithFixup(ntpManipulator, thisVtxPlane, thisVtxZ);

  double thisEvtTime; double thisToFEventTime;
  GetStripEventTime(ntpEvent, thisVtxPlane, thisVtxZ,
                    thisEvtTime, thisToFEventTime);

  // Smallest difference in event time to another event
  Double_t timeDiff = ANtpDefVal::kFloat;
  // Index of the event that had the lowest timeDiff
  Int_t closeEventIndex = ANtpDefVal::kInt;
  // Smallest different in vtx z to another event
  Float_t minDeltaZ = ANtpDefVal::kFloat;
  // Minimum time-of-flight-corrected time difference
  Double_t mintime = ANtpDefVal::kFloat;
 
  // Loop over all the other events to find the nearest in time and space
  for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
    // don't want to compare time to itself
    if(i == event) continue;

    evtManip->SetEventInSnarl(i);
    const NtpSREvent* otherEvent = evtManip->GetEvent();

    int otherVtxPlane; float otherVtxZ;
    GetEvtVtxWithFixup(ntpManipulator, otherVtxPlane, otherVtxZ);

    double otherEvtTime; double otherToFEventTime;
    GetStripEventTime(otherEvent, otherVtxPlane, otherVtxZ,
                      otherEvtTime, otherToFEventTime);
    // time of flight correction for second event
    if( otherToFEventTime !=ANtpDefVal::kDouble) 
      otherToFEventTime -= ((otherVtxZ-thisVtxZ)/Munits::c_light);

    const double deltaT = otherEvtTime - thisEvtTime;
    if(TMath::Abs(deltaT) < TMath::Abs(timeDiff)){
      timeDiff = deltaT;
      closeEventIndex = i;
      eventInfo->closeTimeDeltaZ = thisVtxZ - otherVtxZ;
    }

    // find minimum DeltaZ
    const float deltaZ = thisVtxZ - otherVtxZ;
    if(TMath::Abs(deltaZ) < TMath::Abs(minDeltaZ)){
      minDeltaZ = deltaZ;
    }

    if(thisToFEventTime != ANtpDefVal::kDouble  &&  
       otherToFEventTime != ANtpDefVal::kDouble)
      if(fabs(thisToFEventTime-otherToFEventTime) < fabs(mintime))
        mintime = thisToFEventTime-otherToFEventTime;
  } // end for i

  eventInfo->minTimeSeparation = timeDiff;
  eventInfo->closeTimeEvent = closeEventIndex;
  eventInfo->minDeltaZ = minDeltaZ;
  eventInfo->medianTime = thisEvtTime;
  eventInfo->evtTimeDiff = mintime;
  
  // ADJACENT STRIP VARIABLES

  // Map from plane number to set of strip numbers in that plane
  map<int, set<int> > mplanes;
  // Fill map of strips for all other events in this snarl
  for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
    // don't want to compare time to itself
    if(i == event) continue;
    evtManip->SetEventInSnarl(i);
    const NtpSREvent* otherEvent = evtManip->GetEvent();

    // Fill mplanes
    for(Int_t ns = 0; ns < otherEvent->nstrip; ++ns){
      const int idx = otherEvent->stp[ns];
      if(idx < 0) continue;
      NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(idx);

      if(ntpStrip->ph1.pe >= 2)
        mplanes[ntpStrip->plane].insert(ntpStrip->strip);
    } // end for ns
  } // end for i

  // Number of strips in the event
  int nstrips =0;
  // Count number of strips with neighbouring strips
  int sharedst=0;

  // loop over strips in current event
  for(int ns = 0; ns < ntpEvent->nstrip; ++ns){
    const int stripIdx = ntpEvent->stp[ns];
    if(stripIdx < 0) continue;
    NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(stripIdx);
 
    if(ntpStrip->ph1.pe < 2.) continue;
    ++nstrips;
    map<int, set<int> >::const_iterator PlanesIter =
      mplanes.find(ntpStrip->plane);
    
    //if no hits on same plane, go to next strip
      if(PlanesIter == mplanes.end()) continue;
      set<int> stps = PlanesIter->second;
      const unsigned int stripn = ntpStrip->strip;
      
      //when I fill histograms i want only min separation per event
      for(set<int>::iterator stpit = stps.begin();
          stpit != stps.end();
          ++stpit){
        // Only neighbouring strips
        if(abs(int(stripn)-int(*stpit)) <= 1){
          ++sharedst;
          break;
        }
      } // end for stpit
    } // end for ns

  if(nstrips) 
    eventInfo->sharedStripFraction = float(sharedst)/nstrips;
  
  // SLICE PH FRACTION

  ntpManipulator->GetSliceManipulator()->SetSliceInSnarl(ntpEvent->slc);
  NtpSRSlice* ntpSlice = ntpManipulator->GetSliceManipulator()->GetSlice();

  float slicePH= 0;
  float eventPH = 0;
  for(int it = 0; it < ntpSlice->nstrip; ++it){
    const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpSlice->stp[it]);
    if(ntpStrip->ph1.pe < 2) continue;
    slicePH += ntpStrip->ph1.sigcor;
  }

  for(int it = 0; it < ntpEvent->nstrip; ++it){
   const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpEvent->stp[it]);
   if(ntpStrip->ph1.pe < 2) continue;
   eventPH += ntpStrip->ph1.sigcor;
  }

  if(slicePH)
    eventInfo->slicePHFraction = eventPH/slicePH;

  //NUMBER OF CONSECUTIVE PLANES
  Int_t consecutive=1;
  vector <int> stripplanes;
  for(int it = 0; it < ntpEvent->nstrip; ++it){
    const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpEvent->stp[it]);
    if(ntpStrip->ph1.pe<2.) continue;
    
    bool flag=false;
    for(unsigned int iv=0;iv<stripplanes.size();iv++)
      if(stripplanes[iv]==ntpStrip->plane) flag=true;
    if(!flag) stripplanes.push_back(ntpStrip->plane);
  }
  
  for(unsigned int iv=0;iv<stripplanes.size();iv++)
    sort(stripplanes.begin(),stripplanes.end());
  
  for(unsigned int iv=1;iv<stripplanes.size();iv++)
    if((-stripplanes[0]+stripplanes[iv]-(consecutive-1))==1)
      consecutive++;
  
  eventInfo->consecutivePlanes=consecutive;
  
  // EDGE ACTIVITY VARIABLES

  const double activityTimeStart = thisEvtTime - 4e-8; // +/- 40ns window
  const double activityTimeStop  = thisEvtTime + 4e-8;

  Int_t edgeActivityStrips=0;
  Float_t edgeActivityPH=0;
  Int_t oppEdgeStrips=0;
  Float_t oppEdgePH=0;

  // loop over strips in snarl
  for(int i = 0; i <= fStripArray->GetLast(); ++i){
    const NtpSRStrip* strip = (NtpSRStrip*)fStripArray->At(i);

    // consider only strips in time window
    if(strip->time1 <= activityTimeStart || strip->time1 >= activityTimeStop)
      continue;

    // Only want calorimeter strips
    if(strip->plane >= 121) continue;

    if(strip->ph1.pe < 2) continue;

      // look at U planes in calorimeter
      if(strip->plane % 2 == 1 && strip->tpos < -0.24){
        ++edgeActivityStrips;
        edgeActivityPH += strip->ph1.sigcor;
      }
      // look at opposite edge in U (use 3 strips of fully instr.)
      if(strip->plane % 2 == 1 && strip->tpos > 2.27){
        ++oppEdgeStrips;
        oppEdgePH += strip->ph1.sigcor;
      }
      // look at V planes in calorimeter
      if(strip->plane % 2 == 0 && strip->tpos > 0.24){
        ++edgeActivityStrips;
        edgeActivityPH += strip->ph1.sigcor;
      }
      // look at opposite edge in V (use 3 strips of fully instr.)
      if(strip->plane % 2 == 0 && strip->tpos <- 2.27){
        ++oppEdgeStrips;
        oppEdgePH += strip->ph1.sigcor;
      }
   } // end for i

  eventInfo->edgeActivityPH = edgeActivityPH;
  eventInfo->edgeActivityStrips = edgeActivityStrips;
  eventInfo->oppEdgePH = oppEdgePH;
  eventInfo->oppEdgeStrips = oppEdgeStrips;

  // LARGEST EVENT

  int largestIndex = 0;
  double largestPH = 0;
  for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
    evtManip->SetEventInSnarl(i);
    const NtpSREvent* otherEvent = evtManip->GetEvent();
    if(otherEvent->ph.sigcor > largestPH){
      largestPH = otherEvent->ph.sigcor;
      largestIndex = i;
    }
  }

  if(eventInfo->index == largestIndex)
    eventInfo->largestEventInSnarl = 1;

  //make sure to set the event manipulator back to the correct event
  ntpManipulator->GetEventManipulator()->SetEventInSnarl(event);
}


//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::InitializekNN(ANtpRecoNtpManipulator *ntpManipulator)
{

  if(!fkNNSet){

    //set the location of the weight files based on the detector
    TString weightFileName = "knn.train.far.cedar.daikon.root";
    if(fDetector == Detector::kNear) weightFileName = "knn.train.near.cedar.daikon.root";

    TString base = getenv("SRT_PRIVATE_CONTEXT");
    TString ncutils = "/NCUtils/data/";
    if(base != "" && base != "."){
      //check that the private context director exists, if not use public context
      void *dirptr = gSystem->OpenDirectory(base+ncutils);
      if(!dirptr){
        base = getenv("SRT_PUBLIC_CONTEXT");
      }//end if private context doesnt exist
    }//end if private context is defined
    else base = getenv("SRT_PUBLIC_CONTEXT");

    //check that the public context exists, if not assert false
    if(base == ""){
      MSG("ANtpInfoObjectFillerNC", Msg::kFatal) << "no SRT_PUBLIC_CONTEXT set"
                                                 << endl;
      assert(false);
    }

    TString weightFilePath = base+ncutils+weightFileName;

    //   MSG("ANtpInfoObjectFillerNC", Msg::kInfo) << "kNN weight files "
    //                                      << weightFilePath << endl;
    TString baseConf = getenv("SRT_PUBLIC_CONTEXT");

    Registry reg;
    reg.Set("InterfaceConfigPath", baseConf+"/PhysicsNtuple/Config/Config2007Real.txt");
    reg.Set("FillkNNFilePath", weightFilePath);

    //These are the most recently trained files to use but only exist at fermilab
    //If you want run at fermilab uncomment the next 3 lines

    //reg.Set("FillkNNFilePath", "/afs/fnal.gov/files/home/room1/rustem/data/muon-knn/knn.physics.near.daikon_00.cedar_phy.L010z185i.root");
    //if(fDetector == Detector::kFar)
    //reg.Set("FillkNNFilePath", "/afs/fnal.gov/files/home/room1/rustem/data/muon-knn/knn.physics.far.daikon_00.cedar_phy.L010z185i.root");

    reg.LockValues();

    fAnpInterface = new Anp::Interface();
    fAnpInterface->Config(reg);
    fkNNSet = true;
  }
  fAnpInterface->FillSnarl(ntpManipulator->GetNtpStRecord());

  return;

}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a track
void ANtpInfoObjectFillerNC::FillTrackInformation(ANtpRecoNtpManipulator *ntpManipulator,
                                                  NtpSRTrack *ntpTrack,
                                                  NtpSREvent *ntpEvent,
                                                  ANtpTrackInfoNC *trackInfo,
                                                  ANtpEventInfoNC *eventInfo)
{

  trackInfo->Reset();
  ANtpInfoObjectFiller::FillTrackInformation(ntpTrack, trackInfo);

  if(trackInfo->totalStrips > 0){
    trackInfo->phPerStrip = trackInfo->pulseHeight/(1.*trackInfo->totalStrips);
    trackInfo->phPerPlane = trackInfo->pulseHeight/(1.*trackInfo->planes);
  }

  NtpSRStrip *ntpStrip = 0;

  Int_t plane = 0, strip = 0, stpCtr = 0;
  Float_t xPos = 0., yPos = 0., zPos = 0.;
  Int_t uXTalkStripCtr = 0;
  Int_t vXTalkStripCtr = 0;

  Int_t numDigits = 0, numDoubleEndedStrips = 0;
  Int_t stripsInPlane[500], trkStripsInPlane[500];

  //....................................................................
  // Estimate dtdz of track, using linear regression of (y = Gx +H)
  // <G> = (<xy> - <x><y>) / (<x^2> - <x>^2)
  Double_t sumZT(0.); Double_t sumZ2(0.);
  Double_t sumT(0.); Double_t sumZ(0.);
  Int_t nHits(0); // # of digits at FD, # of strips at ND. whew!

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

    ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrackNS();

    trackInfo->aTrkpass_ns         =ntpTrack->fit.pass;
    trackInfo->aTrkph_ns           =ntpTrack->ph.sigcor;
    trackInfo->aTrklen_ns          =ntpTrack->plane.ntrklike;
    if(ntpTrack->plane.ntrklike>0) trackInfo->aTrkphperpl_ns   = ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
    if(ntpEvent->ph.sigcor>0)      trackInfo->aTrkphper_ns     = ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
    trackInfo->aTrkplu_ns          =ntpTrack->plane.nu;
    trackInfo->aTrkplv_ns          =ntpTrack->plane.nv;
    trackInfo->aTrkstp_ns          =ntpTrack->nstrip;
    trackInfo->aTrkvtx_ns          =ntpTrack->vtx.z;

  for(Int_t ii = 0; ii < 500; ++ii){
    stripsInPlane[ii] = 0;
    trkStripsInPlane[ii] = 0;
  }

  //loop over the strips in the record and keep track of how many
  //are in each plane
  Float_t totalPH = 0.;

  for(Int_t i = 0; i < ntpEvent->nstrip; ++i){

    if(ntpEvent->stp[i] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[i]));
    else continue;
    plane = (int)ntpStrip->plane;
    strip = ntpStrip->strip;
    totalPH = 0;

    if( fDetector == Detector::kFar
        && !fStripIsXTalkEast[plane][strip] ) totalPH += ntpStrip->ph0.sigcor;
    if( !fStripIsXTalkWest[plane][strip] ) totalPH += ntpStrip->ph1.sigcor;

    if(totalPH >= 90.){
      MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
        << "plane " << plane << " strip " << ntpStrip->strip << " is ok "
        << fStripIsXTalkEast[plane][strip] << "/"
        << fStripIsXTalkWest[plane][strip] << " "
        << ntpStrip->ph0.sigcor << "/" << ntpStrip->ph1.sigcor << endl;
      ++stripsInPlane[(int)ntpStrip->plane];
    }

  }//end loop over event strips to see how many are in each plane


  // Now loop over all strips in the track
  // Count the number of strips in each plane for the track
  // Also sum up variables needed to fit t(z)

  for(Int_t i = 0; i < ntpTrack->nstrip; ++i){

    //get the NtpSRStrip object
    if (ntpTrack->stp[i] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpTrack->stp[i]));
    else continue;
    plane = (int)ntpStrip->plane;
    totalPH = 0;

    if( fDetector > 0 && !fStripIsXTalkEast[plane][strip] ) totalPH += ntpStrip->ph0.sigcor;
    if( !fStripIsXTalkWest[plane][strip] ) totalPH += ntpStrip->ph1.sigcor;

    if(totalPH >= 90.)  ++trkStripsInPlane[(int)ntpStrip->plane];

    if(fDetector == Detector::kFar
       && fStripIsXTalkEast[plane][strip] && fStripIsXTalkWest[plane][strip]){
      if((int)ntpStrip->planeview == kU) ++uXTalkStripCtr;
      else if((int)ntpStrip->planeview == kV) ++vXTalkStripCtr;
    }
    else if(fDetector == Detector::kNear && fStripIsXTalkWest[plane][strip]){
      if((int)ntpStrip->planeview == kU) ++uXTalkStripCtr;
      else if((int)ntpStrip->planeview == kV) ++vXTalkStripCtr;
    }
    //..................................................................
    if (ntpStrip->ph0.pe > 2.) {
      sumZ += ntpStrip->z;
      sumZ2 += ntpStrip->z * ntpStrip->z;
      sumZT += ntpStrip->z * ntpStrip->time0;
      sumT += ntpStrip->time0;
      ++nHits;
    }

    if (ntpStrip->ph1.pe > 2.) {
      sumZ += ntpStrip->z;
      sumZ2 += ntpStrip->z * ntpStrip->z;
      sumZT += ntpStrip->z * ntpStrip->time1;
      sumT += ntpStrip->time1;
      ++nHits;
    }
    //..................................................................

  } //end loop over strips in track

  trackInfo->xTalkStrips = uXTalkStripCtr+vXTalkStripCtr;

  //....................................................................
  // Calculate dtdz
  trackInfo->dtdz = 0;
  if (nHits > 1){  // VarZ == 0 if nHits == 1; next conditon would fail
    Double_t CovZT = sumZT - (sumZ*sumT)/(Float_t)nHits;
    Double_t VarZ = sumZ2 - (sumZ*sumZ)/(Float_t)nHits;
    if (VarZ > 0) trackInfo->dtdz = Munits::c_light * CovZT/VarZ;
  }
  //....................................................................


  //loop over all the planes and if the number of strips in the plane for
  //the track is the same as the number for the record, then it is a tracklike plane
  trackInfo->trackLikePlanes = 0;
  for(Int_t i = 0; i < 500; ++i){

    if(stripsInPlane[i]>0 && trkStripsInPlane[i]>0
       && stripsInPlane[i]*0.9 <= trkStripsInPlane[i])
      ++trackInfo->trackLikePlanes;

  }


  Int_t index = 0;
  Int_t nVPlanes = 0, nUPlanes = 0;

  nVPlanes = ntpTrack->plane.nv;
  nUPlanes = ntpTrack->plane.nu;

  trackInfo->uvAsymmetry = 0.;
  if(nVPlanes+nUPlanes>0)
    trackInfo->uvAsymmetry = TMath::Abs((1.*nUPlanes-1.*nVPlanes)/(1.*nVPlanes+1.*nUPlanes));

  trackInfo->signalUseFraction = 0.;
  if(eventInfo->pulseHeight>0.)
    trackInfo->signalUseFraction = trackInfo->pulseHeight/eventInfo->pulseHeight;

  trackInfo->planeUseFraction = 0.;
  if(eventInfo->planes>0)
    trackInfo->planeUseFraction = 1.*trackInfo->trackLikePlanes/(1.*eventInfo->planes);

  trackInfo->totalStrips = 0;

  //does the pathlength of the track make sense?
  if(trackInfo->length<50.){

    //loop over all strips in the track
    trackInfo->totalStrips = (int)ntpTrack->nstrip;
    numDoubleEndedStrips = 0;

    stpCtr = 0;
    //ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(ntpTrack->stp[0]));
    for(Int_t ns = 0; ns < trackInfo->totalStrips; ++ns){

      //get the index for this strip
      if (ntpTrack->stp[ns] >= 0) index = ntpTrack->stp[ns];
      else continue;
      //get the NtpSRStrip object
      ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(index));
      plane = (int)ntpStrip->plane;
      strip = (int)ntpStrip->strip;
      numDigits = (int)ntpStrip->ndigit;

      xPos = ntpTrack->stpx[ns];
      yPos = ntpTrack->stpy[ns];
      zPos = ntpTrack->stpz[ns];

      if(numDigits == 2){

        ++numDoubleEndedStrips;

      } // end if the strip had double sided readout
    } // end loop over strips in the track

    trackInfo->twoEndStripFraction = 0.;
    if(trackInfo->totalStrips>0)
      trackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*trackInfo->totalStrips);

  }//end if reasonable track length

  //fill rustem's pid parameter here. Using new interface.
  trackInfo->kNN = fAnpInterface->GetVar("knn_pid",ntpEvent);
  trackInfo->numScintPlanes = fAnpInterface->GetVar("knn_01",ntpEvent);
  trackInfo->meanSigCor = fAnpInterface->GetVar("knn_10",ntpEvent);
  trackInfo->meanLowStripDivHighStrip = fAnpInterface->GetVar("knn_20",ntpEvent);
  trackInfo->trackSigCorFraction = fAnpInterface->GetVar("knn_40",ntpEvent);

  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a shower
void ANtpInfoObjectFillerNC::FillShowerInformation(NtpSRShower *ntpShower,
                                                   ANtpShowerInfoNC *showerInfo,
                                                   ANtpEventInfoNC *eventInfo)
{

  showerInfo->Reset();
  ANtpInfoObjectFiller::FillShowerInformation(ntpShower, showerInfo);

  if(showerInfo->totalStrips > 0){
    showerInfo->phPerStrip = showerInfo->pulseHeight/(1.*showerInfo->totalStrips);
    showerInfo->phPerPlane = showerInfo->pulseHeight/(1.*showerInfo->planes);
  }

  NtpSRStrip *ntpStrip = 0;

  showerInfo->energyGeV = ntpShower->ph.gev;

  Int_t index = 0;
  Int_t numDigits = 0;
  Int_t numDoubleEndedStrips = 0;

  showerInfo->planeUseFraction = -1.;
  showerInfo->signalUseFraction = -1.;
  if(eventInfo->pulseHeight>0.)
    showerInfo->signalUseFraction = showerInfo->pulseHeight/(eventInfo->pulseHeight);
  if(eventInfo->planes>0)
    showerInfo->planeUseFraction = 1.*showerInfo->planes/(1.*eventInfo->planes);

  Int_t plane = 0;
  Int_t strip = 0;

  showerInfo->xTalkStrips = 0;
  showerInfo->twoEndStripFraction = -1.;

  Float_t sumTposU = 0.;
  Float_t sum2TposU = 0.;
  Float_t sumPHU = 0.;
  Float_t maxTposU = -10.;
  Float_t minTposU = 10.;

  Float_t sumTposV = 0.;
  Float_t sum2TposV = 0.;
  Float_t sumPHV = 0.;
  Float_t maxTposV = -10.;
  Float_t minTposV = 10.;

  Float_t sumPHsigcor(0.);
  Int_t nstrip = ntpShower->nstrip;

  MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "shower strips = "
                                             << ntpShower->nstrip << endl;

   showerInfo->aShwdig_ns          = ntpShower->ndigit;
   showerInfo->aShwplu_ns          = ntpShower->plane.nu;
   showerInfo->aShwplv_ns          = ntpShower->plane.nv;
//    showerInfo->aShwph_ns        = ntpShower->ph.sigcor;
//    showerInfo->aShwphper_ns        = ntpShower->ph.sigcor/ntpEvent ->ph.sigcor;
//    showerInfo->aShwphperpl_ns      = ntpShower->ph.sigcor/ntpShower->plane.n;
//    showerInfo->aShwphperdig_ns     = ntpShower->ph.sigcor/ntpShower->ndigit;
//    showerInfo->aShwphperstp_ns     = ntpShower->ph.sigcor/ntpShower->nstrip;


  for(Int_t ns = 0; ns < ntpShower->nstrip; ++ns){

    //get the index for this strip
    index = ntpShower->stp[ns];

    //get the NtpSRStrip object
    ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(index));
    plane = (int)ntpStrip->plane;
    strip = ntpStrip->strip;

    //......................................................................
    // -=NEW VARIABLES=-
    // use all strips for now
    sumPHsigcor += (ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
    //......................................................................

    numDigits = (int)ntpStrip->ndigit;
    if(numDigits == 2) ++numDoubleEndedStrips;

    if(fDetector == Detector::kNear){
      if(fStripIsXTalkWest[plane][strip]) ++showerInfo->xTalkStrips;
    }
    else if(fDetector == Detector::kFar){
      if(fStripIsXTalkEast[plane][strip] && fStripIsXTalkWest[plane][strip])
        ++showerInfo->xTalkStrips;
    }//end if far detector

    if(ntpStrip->ndigit == 1) continue;
    Float_t stpPH(ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
    if (ntpStrip->planeview == PlaneView::kU){
      maxTposU = TMath::Max(maxTposU, ntpStrip->tpos);
      minTposU = TMath::Min(minTposU, ntpStrip->tpos);
      sumTposU += (ntpStrip->tpos)*stpPH;
      sum2TposU += (ntpStrip->tpos * ntpStrip->tpos) * stpPH;
      sumPHU += stpPH;
    }
    else if (ntpStrip->planeview == PlaneView::kV){
      maxTposV = TMath::Max(maxTposV, ntpStrip->tpos);
      minTposV = TMath::Min(minTposV, ntpStrip->tpos);
      sumTposV += ntpStrip->tpos * stpPH;
      sum2TposV += (ntpStrip->tpos * ntpStrip->tpos) * stpPH;
      sumPHV +=stpPH;
    }

  }//end loop over strips

  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "total strips = " << showerInfo->totalStrips << endl;

  if(showerInfo->totalStrips>0)
    showerInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*showerInfo->totalStrips);

  if(sumPHU>0){
    Float_t centroidU = sumTposU;
    if(sumPHU > 0.) centroidU /= sumPHU;
    else centroidU = 0.;
    showerInfo->transverseRMSU = TMath::Sqrt(TMath::Abs(sum2TposU/sumPHU - centroidU*centroidU));
  }
  else showerInfo->transverseRMSU = 0.;

  if(sumPHV>0){
    Float_t centroidV = sumTposV;
    if(sumPHV > 0.) centroidV /= sumPHV;
    else centroidV = 0.;
    showerInfo->transverseRMSV = TMath::Sqrt(TMath::Abs(sum2TposV/sumPHV - centroidV*centroidV));
  }
  else showerInfo->transverseRMSV = 0.;

  //......................................................................
  // -=NEW VARIABLES=-
  // fills the variables emFrac and phKurtosis

    Float_t meanPHsigcor = sumPHsigcor/nstrip;
    Float_t secondMoment(0.);
    Float_t fourthMoment(0.);

    for (Int_t stp_index_shw = 0; stp_index_shw < nstrip ; ++stp_index_shw){
      Int_t stp_index = ntpShower->stp[stp_index_shw];
      NtpSRStrip* strip =
        dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_index));
      Float_t stpPH(strip->ph0.sigcor + strip->ph1.sigcor);
      secondMoment += TMath::Power((stpPH - meanPHsigcor),2);
      fourthMoment += TMath::Power((stpPH - meanPHsigcor),4);
    }
    Float_t sigma = secondMoment/(nstrip-1);

    // Kurtosis of ph distribution as defined in Numerical Recipes
    showerInfo->phKurtosis = fourthMoment/(nstrip*sigma*sigma)-3.0;

    Int_t emLikeStrips(0);
    Bool_t cluFailure(false); // check for very rare cluster failure
    for (Int_t clu_index_shw = 0; clu_index_shw < ntpShower->ncluster ;
         ++clu_index_shw){
      Int_t clu_index = ntpShower->clu[clu_index_shw];

      if (clu_index<0) {
        MSG("ANtpInfoObjectFillerNC",Msg::kError)
          << "Array index for cluster less than 0. This should never happen!"
          << endl;
        cluFailure = true;
        continue;
      }

      NtpSRCluster* cluster =
        dynamic_cast<NtpSRCluster *>(fClusterArray->At(clu_index));
      if (cluster->id == 0)
        emLikeStrips += cluster->nstrip;
    }
    if (!cluFailure)
      showerInfo->emFrac = (1.0* emLikeStrips)/ntpShower->nstrip;
    else
      showerInfo->emFrac = -1;

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

  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a showers and tracks
void ANtpInfoObjectFillerNC::FillCrossOverInformation(NtpSRTrack *ntpTrack,
                                                      NtpSRShower *ntpShower,
                                                      NtpSREvent *ntpEvent,
                                                      ANtpRecoNtpManipulator *ntpManipulator,
                                                      ANtpTrackInfoNC *trackInfo,
                                                      ANtpShowerInfoNC *showerInfo,
                                                      ANtpEventInfoNC *eventInfo)
{

  Int_t indexShowerStrip = 0;
  Int_t indexTrackStrip = 0;

  eventInfo->trackStripsInShower = 0;
  eventInfo->trackExtension = 0;

  if(ntpShower && ntpTrack && trackInfo->length < 50.){
    for(Int_t ns = 0; ns < ntpShower->nstrip; ++ns){

      //get the index for this strip
      indexShowerStrip = ntpShower->stp[ns];

      //loop over the track strips
      for(Int_t nt = 0; nt < ntpTrack->nstrip; ++nt){

        indexTrackStrip = ntpTrack->stp[nt];

        //check to see if the shower index matches the track index, if so,
        //increment the counter for the number of track strips in the shower
        if(indexTrackStrip == indexShowerStrip) ++eventInfo->trackStripsInShower;

      }//end loop over track strips

    }//end loop over shower strips

    if(ANtpDefaultValue::IsDefault(showerInfo->planes) ||
       ANtpDefaultValue::IsDefault(trackInfo->planes) )
      eventInfo->trackExtension = ANtpDefaultValue::kInt;
    else
      eventInfo->trackExtension = showerInfo->planes - trackInfo->planes;
  }//end if shower and track are both present and track length < 50. m

  //fill niki's variables
  Int_t    ssind,ssplind;
  Double_t ssphtot;
  Bool_t   foundshw,foundtrk;

  Int_t planes = ntpEvent->plane.beg;

  Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
  ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;

  NtpSRTrack  *ntpTrackNS = ntpManipulator->GetEventManipulator()->GetPrimaryTrackNS();

  // loop over strips to compute ph3 ph6 phlast phcommon
  for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
    Int_t stp_index=((ntpEvent->stp)[evss]);

    if(stp_index!=-1){

     NtpSRStrip *ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_index));
     if(!ntpStrip) continue;
     ssind   = ntpStrip->strip;
     ssplind = ntpStrip->plane;
     ssphtot = ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;

     foundshw=false;
     foundtrk=false;

     Double_t phstrips=0;
     Double_t phstript=0;
    // shower strips
     Int_t sshwind,sshwplind;
     Double_t sshwphtot;

    if(ntpEvent->nshower > 0) {
      for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
       Int_t stp_indexs=((ntpShower->stp)[jj]);

       if(stp_indexs!=-1){

        NtpSRStrip *ntpStrip1 = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_indexs));
        if(!foundshw){
          sshwind   =ntpStrip1->strip;
          sshwplind =ntpStrip1->plane;
          sshwphtot =ntpStrip1->ph1.sigcor+ntpStrip1->ph0.sigcor;
          if(sshwind==ssind && sshwplind==ssplind) {
            foundshw=true;
            phstrips=sshwphtot;
          }
        }
       }
      }
    }
    // tracks strips
    Int_t strkind,strkplind;
    Double_t strkphtot;

    if(ntpEvent->ntrack > 0) {

      for(Int_t jj=0;jj<ntpTrackNS->nstrip;jj++){
        Int_t stp_indext=((ntpTrackNS->stp)[jj]);

        if(stp_indext!=-1){

        NtpSRStrip *ntpStript = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_indext));

        if(!foundtrk){
          strkind  =ntpStript->strip;
          strkplind=ntpStript->plane;
          strkphtot=ntpStript->ph1.sigcor+ntpStript->ph0.sigcor;
          if(strkind==ssind && strkplind==ssplind) {
            foundtrk=true;
            phstript=strkphtot;
          }
        }

       }
      }
    }

    if(foundshw && foundtrk) {
      hitcommon= hitcommon+1;
      phcommon = phcommon+phstrips+phstript;
    }
    if(!foundshw && ! foundtrk) {
      hitnowere=hitnowere+1;
      phnowere=phnowere+ssphtot;
    }
    if(ssplind>=planes && ssplind<=(planes+3)){
      ph3=ph3+ssphtot;
    }
    else if(ssplind>(planes+3) && ssplind<=(planes+6)){
      ph6=ph6+ssphtot;
    }
    else {
      phlast=phlast+ssphtot;
    }

   } // end of if stp_index != -1

  } // end of looping over event strips....

//
   eventInfo->aPh3_ns       = ph3;
   eventInfo->aPh6_ns       = ph6;
   eventInfo->aPhlast_ns    = phlast;
   eventInfo->aPhcommon_ns  = phcommon;


  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to an event
void ANtpInfoObjectFillerNC::FillEventInformation(ANtpRecoNtpManipulator *ntpManipulator,
                                                  NtpSREvent *ntpEvent,
                                                  ANtpEventInfoNC *eventInfo)
{

  eventInfo->Reset();
  ANtpInfoObjectFiller::FillEventInformation(ntpManipulator,
                                             ntpEvent,
                                             eventInfo);

  NtpSRStrip *ntpStrip = 0;

  eventInfo->lengthInPlanes = eventInfo->endPlane - eventInfo->begPlane + 1;
  eventInfo->lateBucketPHFraction = ntpEvent->bleach.lateBucketPHFraction;
  eventInfo->timeWeightedPHFraction = ntpEvent->bleach.timeWeightedPHFraction;
  eventInfo->straightPHFraction = ntpEvent->bleach.straightPHFraction;
  eventInfo->fixedWindowPH = ntpEvent->bleach.fixedWindowPH;

  //add up the pulse height in each plane to find the maximum in the event
  Float_t planePH[500] = {0.};
  eventInfo->totalStrips = ntpEvent->nstrip;
  eventInfo->passStrips = 0;
  eventInfo->modifiedPH = 0.;
  if(eventInfo->totalStrips > 0){
    eventInfo->phPerStrip = eventInfo->pulseHeight/(1.*eventInfo->totalStrips);
    eventInfo->phPerPlane = eventInfo->pulseHeight/(1.*eventInfo->planes);
  }

  double maxStripTime = -1.e12;
  double minStripTime = 1.e12;

//......................................................................
  // -=NEW VARIABLES=-

  // PLANE ADC (`GLOBAL') VARIABLES
  // First make summary of plane ADC
  const UInt_t nPlanes(eventInfo->planes);
  vector<Float_t> planeAdc(nPlanes+12,0);
  // pad vector with 12 zeros. these are needed later
  const vector<Float_t>::const_iterator pAdc_begin(planeAdc.begin());
  const vector<Float_t>::const_iterator pAdc_last(pAdc_begin + nPlanes);

  // Also calculate some simple event duration type variables
  Double_t t_1stStp(10.); //time of first strip in event
  Double_t t_lastStp(0.); //  & of last
  Double_t t_sumStp(0.);  // sum of times
  Double_t t_sum2Stp(0.); // sum of squares of times
  Int_t nStpHit( ntpEvent->nstrip ); // number of double ended strips
//......................................................................

  for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){

    //get the current strip
    if (ntpEvent->stp[ns] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[ns]));
    else continue;
    if(ntpStrip->plane>0 && ntpStrip->plane<486){
      planePH[ntpStrip->plane] += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;

      if( fDetector == Detector::kFar
          && (!fStripIsXTalkWest[ntpStrip->plane][ntpStrip->strip]
              || !fStripIsXTalkEast[ntpStrip->plane][ntpStrip->strip]))
        ++eventInfo->passStrips;
      else if( fDetector == Detector::kNear
               && !fStripIsXTalkWest[ntpStrip->plane][ntpStrip->strip])
        ++eventInfo->passStrips;

      if(ntpStrip->ph1.raw > 200.) eventInfo->modifiedPH += ntpStrip->ph1.raw;

      if(ntpStrip->time1 > maxStripTime) maxStripTime = ntpStrip->time1;
      if(ntpStrip->time1 < minStripTime) minStripTime = ntpStrip->time1;

    }
    //......................................................................
    // -=NEW VARIABLES=-

    // Plane ADC vars
    Int_t pln_index = ntpStrip->plane - eventInfo->begPlane;
    if (pln_index < (Int_t) nPlanes && pln_index >= 0) {
      planeAdc[pln_index] += ntpStrip->ph0.sigcor;
      planeAdc[pln_index] += ntpStrip->ph1.sigcor;
    }

    // event duration variables
    Double_t t_stpHitA = TMath::Max(ntpStrip->time0 , ntpStrip->time1);
    Double_t t_stpHitB = TMath::Min(ntpStrip->time0 , ntpStrip->time1);
    //Earlier hit may be -ive, meaning a single ended hit
    t_1stStp = TMath::Min(t_1stStp, t_stpHitA);
    t_lastStp = TMath::Max(t_lastStp, t_stpHitA);
    t_sumStp += t_stpHitA;
    t_sum2Stp += t_stpHitA * t_stpHitA;

    if (t_stpHitB > 0 ){
      t_1stStp = TMath::Min(t_1stStp, t_stpHitB);
      t_lastStp = TMath::Max(t_lastStp, t_stpHitB);
      t_sumStp += t_stpHitB;
      t_sum2Stp += t_stpHitB * t_stpHitB;
      ++nStpHit;
    }
    //......................................................................

  }//end loop over event strips

  //......................................................................
  // -=NEW VARIABLES=-
  eventInfo->stripTimeMean = t_sumStp / nStpHit;
  Float_t t_variance =  t_sum2Stp / nStpHit
    - (t_sumStp / nStpHit) * (t_sumStp / nStpHit);

  eventInfo->stripTimeRMS = TMath::Sqrt(TMath::Abs(t_variance));

  eventInfo->stripTime1st = t_1stStp;
  eventInfo->stripTimelast = t_lastStp;
  //......................................................................

  eventInfo->eventDuration = maxStripTime - minStripTime;

  if(fDetector == Detector::kNear) FindEarlyActivityWeight(ntpEvent, eventInfo);

  //find the rms of the ph per plane for the event and the
  //plane with the max ph.  also find the maximum signal in windows of
  //3,6,9,12 planes from the first plane

  TH1F planePHDist("planePHDist", "", 400, 0., 400.);
  eventInfo->maxPlanePH = 0.;
  eventInfo->maxPHIn3Planes = 0.;
  eventInfo->maxPHIn6Planes = 0.;
  eventInfo->maxPHIn9Planes = 0.;
  eventInfo->maxPHIn12Planes = 0.;

  for(Int_t i = eventInfo->begPlane; i <= eventInfo->endPlane; ++i){
    if(planePH[i] > eventInfo->maxPlanePH) eventInfo->maxPlanePH = planePH[i];
    if(planePH[i] > 0.) planePHDist.Fill(planePH[i]*1.e-4);
  }

  Float_t phInWindow = 0.;
  for(Int_t windowSize = 3; windowSize<14; windowSize +=3){
    for(Int_t i = eventInfo->begPlane; i<=eventInfo->endPlane; ++i){
      phInWindow = 0.;
      for(Int_t j = 0; j<windowSize && j+i<=eventInfo->endPlane; ++j){
        phInWindow += planePH[i+j];
      }
      if(windowSize == 3 && eventInfo->maxPHIn3Planes<phInWindow)
        eventInfo->maxPHIn3Planes = phInWindow;
      else if(windowSize == 6 && eventInfo->maxPHIn6Planes<phInWindow)
        eventInfo->maxPHIn6Planes = phInWindow;
      else if(windowSize == 9 && eventInfo->maxPHIn9Planes<phInWindow)
        eventInfo->maxPHIn9Planes = phInWindow;
      else if(windowSize == 12 && eventInfo->maxPHIn12Planes<phInWindow)
        eventInfo->maxPHIn12Planes = phInWindow;
    } //end loop over planes in event

  }//end loop over window size

  eventInfo->pulseHeightRms = planePHDist.GetRMS();

  //......................................................................
  // -=NEW VARIABLES=-

  //Calculate pulse height in plane groups 1 to 3 and 4 to 6
  eventInfo->triPlane1PH = accumulate(pAdc_begin, pAdc_begin+3, 0.0 );
  eventInfo->triPlane2PH = accumulate(pAdc_begin+3, pAdc_begin+6, 0.0 );
  if (nPlanes>6)
    eventInfo->triPlaneOverPH = accumulate(pAdc_begin+6, pAdc_last, 0.0 );
  else
    eventInfo->triPlaneOverPH = 0.0;
  //......................................................................

  eventInfo->eventSummaryPlanes = ntpManipulator->GetSnarlEventSummary().plane.n;

  // Image recognition discriminants (KEA)
  //FillVHSEvtInfo(ntpManipulator, eventInfo);

  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to the MC truth
void ANtpInfoObjectFillerNC::FillMCInformation(NtpMCTruth *ntpMCTruth,
                                               TClonesArray *stdHepArray,
                                               ANtpTruthInfoBeam *truthInfo)
{
  truthInfo->Reset();
  ANtpInfoObjectFiller::FillMCTruthInformation(ntpMCTruth, truthInfo);
  ANtpInfoObjectFillerBeam::FillBeamMCTruthInformation(ntpMCTruth, stdHepArray, truthInfo);

  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::FillPlanePixelSignalArrays(NtpSREvent *ntpEvent)
{
  MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "in FillPlanePixelSignalArrays" << endl;

  NtpSRStrip *ntpStrip = 0;

  Int_t pixel = 0;
  Int_t plane = 0;

  //fill the plane to pmt map first.  the idea is that there are at most 2 PMT's per plane
  //in either detector.  put one NtpSRStrip->pmtindex0(1) into the fPlaneToPMTMapEast(West)
  //array.  near detector only has pmt's on the west side.  call that pmt 0 for the plane.
  //then when filling the signal arrays if the pmtindex0(1) for a strip is not the one in the
  //array, it is pmt 1 for the plane.
  for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){

    //get the NtpSRStrip object
    if (ntpEvent->stp[ns] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[ns]));
    else continue;
    plane = ntpStrip->plane;

    MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
      << "plane " << plane << " strip " << ntpStrip->strip
      << " pmt1 " << ntpStrip->pmtindex1 << endl;
    MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
      << " pmt0 " << ntpStrip->pmtindex0 << endl;

    fPlaneToPMTMapWest[plane] = ntpStrip->pmtindex1;
    if(fDetector == Detector::kFar) fPlaneToPMTMapEast[plane] = ntpStrip->pmtindex0;

  }//end loop over strips to fill plane to pmt maps

  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "finish filling plane to pmt maps" << endl;

  //loop over all the strips and put the signal size in the arrays
  for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){

    //get the NtpSRStrip object
    if(ntpEvent->stp[ns] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[ns]));
    else continue;
    plane = ntpStrip->plane;

    MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
      << "plane " << plane << " strip " << ntpStrip->strip
      << " detector " << fDetector << endl;

    //got the strip, now fill the plane pixel info
    if(fDetector == Detector::kFar){

      //east first
      pixel = fStripToPixelMapEast[ntpStrip->strip];

      if(ntpStrip->pmtindex0 == fPlaneToPMTMapEast[plane])
        fPlanePixelEastSignal[plane][pixel][0] = ntpStrip->ph0.sigcor;
      else
        fPlanePixelEastSignal[plane][pixel][1] = ntpStrip->ph0.sigcor;

      //now west
      pixel = fStripToPixelMapWest[ntpStrip->strip];
      if(ntpStrip->pmtindex1 == fPlaneToPMTMapWest[plane])
        fPlanePixelWestSignal[plane][pixel][0] = ntpStrip->ph1.sigcor;
      else
        fPlanePixelWestSignal[plane][pixel][1] = ntpStrip->ph1.sigcor;

    }//end if far detector;
    else if(fDetector == Detector::kNear){

      //figure out what type of plane you are
      //planeview = 2 --> U, 3 --> V
      //map types 2 and 3 are for full planes, 1 is for partial
      MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
        << "in near detector algorithm "
        << " plane " << plane
        << " coverage = " << fPlaneCoverage[plane]
        << " view = " << (int)ntpStrip->planeview << endl;

      if(fPlaneCoverage[plane] == kNearFull){

        if((int)ntpStrip->planeview == kU){
          pixel = fStripToPixelMapNearU2[ntpStrip->strip];
          if(pixel < 0) pixel = fStripToPixelMapNearU3[ntpStrip->strip];
          if(pixel < 0)
            MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
              << "cannot find pixel for strip " << ntpStrip->strip
              << " plane " << plane << endl;

        }
        else if((int)ntpStrip->planeview == kV){
          pixel = fStripToPixelMapNearV2[ntpStrip->strip];
          if(pixel < 0) pixel = fStripToPixelMapNearV3[ntpStrip->strip];
          if(pixel < 0)
            MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
              << "cannot find pixel for strip " << ntpStrip->strip
              << " plane " << plane << endl;

        }

      }//end if full
      else if(fPlaneCoverage[plane] == kNearPartial){

        if((int)ntpStrip->planeview == kU)
          pixel = fStripToPixelMapNearU1[ntpStrip->strip];
        else if((int)ntpStrip->planeview == kV)
          pixel = fStripToPixelMapNearV1[ntpStrip->strip];

      }//end if partial

      if(ntpStrip->pmtindex1 == fPlaneToPMTMapWest[plane])
        fPlanePixelWestSignal[plane][pixel][0] = ntpStrip->ph1.sigcor;
      else
        fPlanePixelWestSignal[plane][pixel][1] = ntpStrip->ph1.sigcor;


    }//end if near detector
  }//end loop over strips

  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::FindEarlyActivityWeight(NtpSREvent *ntpEvent,
                                                     ANtpEventInfoNC *eventInfo)
{

  //find the time of the earliest strip in the event
  double earliestEventTime = 1.e20;
  map<int,int> eventPlanes;
  NtpSRStrip *strip = 0;

  for(int i = 0; i < ntpEvent->nstrip; ++i){
    if (ntpEvent->stp[i] >= 0)
      strip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[i]));
    else continue;
    if(strip->time1 < earliestEventTime) earliestEventTime = strip->time1;
    if(eventPlanes.find(strip->pmtindex1) == eventPlanes.end()) eventPlanes[strip->pmtindex1] = 1;

  }

//   for(int i = 0; i < ntpEvent->nstrip; ++i){
//     strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(ntpEvent->stp[i]));

//     MSG("ANtpInfoObjectFillerNC", Msg::kInfo)
//       << fHeaderInfo->snarl << " " << fHeaderInfo->events
//       << " " << eventInfo->event << " "
//       <<strip->plane << " " << earliestEventTime
//       << " " << 1.e9*(strip->time1-earliestEventTime) << endl;
//   }

  //now find the weigthed sum of ADC for the early strips - make sure the early
  //strips come from the same pmts as the event strips.

  double weightSum = 0.;
  for(int i = 0; i < fStripArray->GetLast()+1; ++i){

    strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(i));

    //work in ns
    if(1.e9*(earliestEventTime - strip->time1) > 0.
       && 1.e9*(earliestEventTime - strip->time1) < 1000.*1.5
       && eventPlanes.find(strip->pmtindex1) != eventPlanes.end()){
      weightSum += strip->ph1.sigcor*TMath::Exp(-1.e9*(earliestEventTime - strip->time1)/700.);

      MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
        << strip->pmtindex1 << " " << strip->plane
        << " " << strip->ph1.sigcor
        << " " << 1.e9*(earliestEventTime-strip->time1)
        << " " << weightSum
        << " " << eventInfo->pulseHeight << endl;

    }
  }

  eventInfo->earlyWeightedADC = weightSum;

  return;
}

//----------------------------------------------------------------------
Float_t ANtpInfoObjectFillerNC::FindNearestNeighborPixelSignal(Int_t plane, Int_t pixel, Int_t pmt,
                                                               Float_t planePixelSignal[][64][2])
{
  MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
    << "in FindNearestNeighborPixelSignal" << endl;

  Float_t sumNeighborSignal = 0.;

  if(fDetector == Detector::kFar){
    if(pixel == 5 || pixel == 6 || pixel == 9 || pixel == 10)
      sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
                           + planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]);
    else if(pixel < 3 && pixel > 0)
      sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
                           + planePixelSignal[plane][pixel+4][pmt]);
    else if(pixel < 15 && pixel > 12)
      sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
                           + planePixelSignal[plane][pixel-4][pmt]);
    else if(pixel > 3 && pixel < 15 && (pixel+1)%4 == 0)
      sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
                           + planePixelSignal[plane][pixel-1][pmt]);

    else if(pixel > 3 && pixel < 12 && pixel%4 == 0)
      sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
                           + planePixelSignal[plane][pixel+1][pmt]);
    else if(pixel == 0)
      sumNeighborSignal = planePixelSignal[plane][pixel+4][pmt] + planePixelSignal[plane][pixel+1][pmt];

    else if(pixel == 3)
      sumNeighborSignal = planePixelSignal[plane][pixel+4][pmt] + planePixelSignal[plane][pixel-1][pmt];

    else if(pixel == 12)
      sumNeighborSignal = planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+1][pmt];

    else if(pixel == 15)
      sumNeighborSignal = planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel-1][pmt];

  }//end if far detector
  else if(fDetector == Detector::kNear){

    if( (pixel > 8 && pixel < 15) || (pixel > 16 && pixel < 23)
        || (pixel > 24 && pixel < 31) || (pixel > 32 && pixel < 39)
        || (pixel > 40 && pixel < 47) || (pixel > 48 && pixel < 55) )
      sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
                           + planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]);
    else if( pixel > 0 && pixel < 7)
      sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
                           + planePixelSignal[plane][pixel+8][pmt]);
    else if( pixel > 56 && pixel < 63)
      sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
                           + planePixelSignal[plane][pixel-8][pmt]);
    else if( pixel > 8 && pixel < 63 && (pixel+1)%8==0 )
      sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
                           + planePixelSignal[plane][pixel-1][pmt]);
    else if( pixel > 7 && pixel < 56 && (pixel)%8==0 )
      sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
                           + planePixelSignal[plane][pixel+1][pmt]);
    else if(pixel==0)
      sumNeighborSignal = planePixelSignal[plane][pixel+8][pmt] + planePixelSignal[plane][pixel+1][pmt];
    else if(pixel==7)
      sumNeighborSignal = planePixelSignal[plane][pixel+8][pmt] + planePixelSignal[plane][pixel-1][pmt];
    else if(pixel==56)
      sumNeighborSignal = planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+1][pmt];
    else if(pixel==63)
      sumNeighborSignal = planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel-1][pmt];
  }

  return sumNeighborSignal;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFillerNC::FindXTalkStrips(NtpSREvent *ntpEvent)
{
  MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "in FindXTalkStrips" << endl;

  NtpSRStrip *ntpStrip = 0;

  Int_t plane = 0;
  Int_t pixelEast = 0;
  Int_t pmtEast = 0;
  Int_t pixelWest = 0;
  Int_t pmtWest = 0;

  for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){

    //get the NtpSRStrip object
    if (ntpEvent->stp[ns] >= 0)
      ntpStrip = dynamic_cast<NtpSRStrip *>
        (fStripArray->At(ntpEvent->stp[ns]));
    else continue;
    plane = (int)ntpStrip->plane;

    if(fDetector == Detector::kFar){
      pixelEast = fStripToPixelMapEast[ntpStrip->strip];
      pmtEast = 0;
      if(fPlaneToPMTMapEast[plane] != ntpStrip->pmtindex0) pmtEast = 1;

      if(ntpStrip->ph0.sigcor < 0.1*FindNearestNeighborPixelSignal(plane, pixelEast,
                                                                   pmtEast, fPlanePixelEastSignal)
         || ntpStrip->ph0.sigcor < 90.){
        fStripIsXTalkEast[plane][ntpStrip->strip] = 1;
        MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
          << "east plane " << plane << " strip " << ntpStrip->strip
          << " is xtalk " << ntpStrip->ph0.sigcor <<  endl;
      }

    }

    pixelWest = fStripToPixelMapWest[ntpStrip->strip];
    pmtWest = 0;
    if(fPlaneToPMTMapWest[plane] != ntpStrip->pmtindex1) pmtWest = 1;

    if(ntpStrip->ph1.sigcor < 0.1*FindNearestNeighborPixelSignal(plane, pixelWest,
                                                                 pmtWest, fPlanePixelWestSignal)
       || ntpStrip->ph1.sigcor < 90.){

      fStripIsXTalkWest[plane][ntpStrip->strip] = 1;
      MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
        << "west plane " << plane << " strip " << ntpStrip->strip
        << " is xtalk " << ntpStrip->ph1.sigcor << endl;
    }

  }//end loop over strips
  return;
}

---