ANtpInfoObjectFiller.cxx

Go to the documentation of this file.

//$Id: ANtpInfoObjectFiller.cxx,v 1.72 2009/11/16 14:40:56 pittam Exp $
//
//ANtpInfoObjectFiller.cxx
//
//B. Rebel 02/2005

#include "AnalysisNtuples/Module/ANtpInfoObjectFiller.h"
#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "AnalysisNtuples/ANtpDefaultValue.h"
#include "MessageService/MsgService.h"
#include "AnalysisNtuples/ANtpEventInfo.h"
#include "AnalysisNtuples/ANtpShowerInfo.h"
#include "AnalysisNtuples/ANtpTrackInfo.h"
#include "AnalysisNtuples/ANtpTruthInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRCosmicRay.h"
#include "CandNtupleSR/NtpSREventSummary.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "BField/BfieldCoilCurrent.h"
#include "DcsUser/CoilTools.h"
#include "DcsUser/Dcs_Mag_Near.h"
#include "DcsUser/BfldDbiCoilState.h"
#include "DcsUser/HvStatusFinder.h"
#include "DcsUser/HvStatus.h"
#include "Conventions/Munits.h"
#include "Conventions/Detector.h"
#include "Conventions/ReleaseType.h"
#include "DataUtil/PlaneOutline.h"
#include "DataUtil/DataQualDB.h"
#include "NtupleUtils/LISieve.h"
#include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
#include "AstroUtil/Ast.h"
#include "AstroUtil/AstTime.h"

#include "TString.h"
#include <cassert>
#include <algorithm>
#include <cmath>

#include <time.h>
static const double kFarSideLength = 1.656854250;    
static const double kInverseSqrt2 = TMath::Power(2., -0.5);

using namespace AstUtil;

ClassImp(ANtpInfoObjectFiller)

CVSID("$Id: ANtpInfoObjectFiller.cxx,v 1.72 2009/11/16 14:40:56 pittam Exp $");

//-------------------------------------------------------------------
ANtpInfoObjectFiller::ANtpInfoObjectFiller() :
  fStripArray(0),
  fDetector(Detector::kUnknown),
  fRowNum(0)
{
    
  MSG("ANtpInfoObjectFiller", Msg::kDebug) 
    << "ANtpInfoObjectFiller::Constructor" << endl;
  
}

//-------------------------------------------------------------------
ANtpInfoObjectFiller::ANtpInfoObjectFiller(Detector::Detector_t detector) :
  fStripArray(0),
  fDetector(detector),
  fRowNum(0)
{
    
  MSG("ANtpInfoObjectFiller", Msg::kDebug) 
    << "ANtpInfoObjectFiller::Constructor" << endl;
  
}

//-------------------------------------------------------------------
ANtpInfoObjectFiller::ANtpInfoObjectFiller(TClonesArray *strips) :
  fStripArray(strips),
  fDetector(Detector::kUnknown),
  fRowNum(0)
{
    
  MSG("ANtpInfoObjectFiller", Msg::kDebug) 
    << "ANtpInfoObjectFiller::Constructor" << endl;
  
}

//----------------------------------------------------------------------
ANtpInfoObjectFiller::~ANtpInfoObjectFiller()
{
  
  MSG("ANtpInfoObjectFiller", Msg::kDebug) 
    << "ANtpInfoObjectFiller::Destructor" << endl;
  
}

//----------------------------------------------------------------------
void ANtpInfoObjectFiller::SetDetector(Detector::Detector_t detector)
{
  fDetector = detector;
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to header information
void ANtpInfoObjectFiller::FillHeaderInformation(ANtpRecoNtpManipulator *ntpManip,
                                                 Detector::Detector_t det,
                                                 SimFlag::SimFlag_t dataType,
                                                 ANtpHeaderInfo *headerInfo)
{
  const NtpSREventSummary& summary(ntpManip->GetSnarlEventSummary());

  headerInfo->Reset();
  
  headerInfo->detector = (int)det;
  headerInfo->dataType = (int)dataType;
  headerInfo->run = ntpManip->GetRun();
  headerInfo->subRun = ntpManip->GetSubRun();
  headerInfo->snarl = ntpManip->GetSnarl();
  //always call it a new snarl because this method is called once per record
  //reset the value in the loop over events in the module
  headerInfo->newSnarl = 1; 
  headerInfo->triggerSource = ntpManip->GetTriggerSource();
  headerInfo->spillType = ntpManip->GetSpillType();
  headerInfo->events = summary.nevent;
  headerInfo->slices = summary.nslice;
  headerInfo->year = (int)summary.date.year;
  headerInfo->month = (int)summary.date.month;
  headerInfo->day = (int)summary.date.day;
  headerInfo->hour = (int)summary.date.hour;
  headerInfo->minute = (int)summary.date.minute;
  headerInfo->utc = summary.date.utc;


  headerInfo->second = summary.date.sec;
  // It sometimes happens in MRCC files that summary.trigtime and
  // summary.date.sec are filled with junk. If this is the case,
  // refill the second based on the utc value, which is still
  // correct. This means we lose the fractional part, but since
  // trigtime is wrong, we couldn't have got it anyway
  if(headerInfo->second<0 || headerInfo->second>61 || isnan(headerInfo->second)){
    MSG("ANtpInfoObjectFiller", Msg::kWarning) 
      << "summary.date.sec out of range (" << summary.date.sec << ") " 
      << "Resetting from utc" << endl;
    // Code shamelessly stolen from VldTimeStamp::GetTime
    time_t atime=summary.date.utc;
    struct tm* ptm=gmtime(&atime);
    headerInfo->second=ptm->tm_sec;
  }


  double longitude = AstUtil::kFarDetLongitude;
  if(det == Detector::kNear) longitude = AstUtil::kNearDetLongitude;
  double gmst = 0.;
  double hour = 1.*headerInfo->hour + 1.*headerInfo->minute/60. + headerInfo->second/3600.;
  AstUtil::CalendarToJulian(headerInfo->year, headerInfo->month, 
                            headerInfo->day, hour, headerInfo->julianDate);
  AstUtil::JulianToGMST(headerInfo->julianDate, gmst);
  AstUtil::GSTToLST(gmst, longitude, headerInfo->localSiderealTime);

  headerInfo->snarlPulseHeight = summary.ph.sigcor;     
  headerInfo->triggerPMTTime = summary.litime;
  // This still won't protect against the case of trigtime being
  // filled with junk that is interpreted as a very tiny positive
  // number, but it's the best we can do I think
  if(summary.trigtime<0 || summary.trigtime>1){
    MSG("ANtpInfoObjectFiller", Msg::kWarning) 
      << "summary.trigtime out of range (" << summary.trigtime << ") " 
      << "Resetting to zero." << endl;
    headerInfo->triggerTime=0;
  }
  else{
    headerInfo->triggerTime = summary.trigtime;
  }
  headerInfo->passedDeMux = 1;
  headerInfo->crateMask = ntpManip->GetDataQuality().cratemask;
  //headerInfo->isGoodFDData = (int)DataUtil::IsGoodFDData(ntpManip->GetNtpStRecord());
  headerInfo->isGoodData = (int)DataUtil::IsGoodData(ntpManip->GetNtpStRecord());

  //......................................................................
  //Now fill new data members
  headerInfo->snarlPE = summary.ph.pe;
  headerInfo->sntpRow = fRowNum++;
  headerInfo->isLIold = LISieve::IsLI( *(ntpManip->GetNtpStRecord()) );
  headerInfo->isLI = LISieve::IsLIforNC( *(ntpManip->GetNtpStRecord()) );

  //......................................................................
  //Fill Release type member
  headerInfo->softVersion = ntpManip->GetReleaseMCType();
  
  //Passed DeMux flag
  if(det == Detector::kFar) {
    const NtpSRDmxStatus& dmx( ntpManip->GetDmxStatusInfo());
    if (dmx.ismultimuon || dmx.validplanesfail) headerInfo->passedDeMux = 0;
  }
  
  //loop over the strips to figure out the 
  double ph2pe = 0.;
  double ph = 0.;
  NtpSRStrip *strip = 0;
  for(int i = 0; i < fStripArray->GetLast()+1; ++i){
    strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(i));
    ph = strip->ph0.pe+strip->ph1.pe;
    if(ph > 2.) ph2pe += strip->ph0.sigcor + strip->ph1.sigcor;
  }
  headerInfo->snarlPulseHeight2PE = ph2pe;
    
  //figure out the coil status
  headerInfo->coilStatus = 0;

  //coil is always on in the MC
  if(dataType == SimFlag::kMC){
    headerInfo->coilStatus = 1;
    return;
  }

  //get the appropriate context
  VldTimeStamp timeStamp(headerInfo->year, headerInfo->month, 
                         headerInfo->day, headerInfo->hour, 
                         headerInfo->minute, 
                         TMath::Nint(headerInfo->second));
  fVldc = VldContext(det, dataType, timeStamp);

  //check out the HV status of the detector - only works for 
  //far detector currently
  if(det == Detector::kFar){
    HvStatus::HvStatus_t hv = HvStatusFinder::Instance().GetHvStatus(fVldc,60,1);
    headerInfo->hvStatus = (int)HvStatus::Good(hv);
    //   MSG("ANtpInfoObjectFiller", Msg::kInfo) << HvStatus::Good(hv)
    //                                    << " " 
    //                                    << headerInfo->hvStatus
    //                                    << endl;
    const BfldDbiCoilState *farstate = CoilTools::Instance().GetCoilState(fVldc);
    if(farstate) headerInfo->coilCurrent =  farstate->GetCurrent();
  }
  else if(det == Detector::kNear){
    const Dcs_Mag_Near *nearstate = CoilTools::Instance().GetMagNear(fVldc);
    if(nearstate) headerInfo->coilCurrent = nearstate->GetCurrent();
  }

  //CoilTools allows for generic access of the database - no 
  //detector dependent stuff needed by user
  if( CoilTools::IsReverse(fVldc) ) headerInfo->coilStatus = -1;
  else headerInfo->coilStatus = 1;
  
  if( !CoilTools::IsOK(fVldc) ) headerInfo->coilStatus = 0;

  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a track
void ANtpInfoObjectFiller::FillTrackInformation(NtpSRTrack *ntpTrack, 
                                                ANtpTrackInfo *trackInfo)
{
  
  trackInfo->Reset();
  trackInfo->planes = ntpTrack->plane.n;
  trackInfo->totalStrips = ntpTrack->nstrip;
  trackInfo->pulseHeight = ntpTrack->ph.sigcor;
  if(TMath::Abs(ntpTrack->momentum.qp) > 0.) 
    trackInfo->fitMomentum = (1./ntpTrack->momentum.qp);
  trackInfo->sigmaQoverP = ntpTrack->momentum.eqp;
  trackInfo->rangeMomentum = ntpTrack->momentum.range;                            
  trackInfo->begPlane = ntpTrack->plane.beg;
  trackInfo->endPlane = ntpTrack->plane.end;
  trackInfo->begPlaneU = ntpTrack->plane.begu;
  trackInfo->endPlaneU = ntpTrack->plane.endu;
  trackInfo->begPlaneV = ntpTrack->plane.begv;
  trackInfo->endPlaneV = ntpTrack->plane.endv;
  trackInfo->length = ntpTrack->ds;
  trackInfo->vtxX = ntpTrack->vtx.x; 
  trackInfo->vtxY = ntpTrack->vtx.y; 
  trackInfo->vtxZ = ntpTrack->vtx.z; 
  trackInfo->dcosXVtx = ntpTrack->vtx.dcosx;
  trackInfo->dcosYVtx = ntpTrack->vtx.dcosy;
  trackInfo->dcosZVtx = ntpTrack->vtx.dcosz;
  trackInfo->endX = ntpTrack->end.x; 
  trackInfo->endY = ntpTrack->end.y; 
  trackInfo->endZ = ntpTrack->end.z; 
  trackInfo->dcosXEnd = ntpTrack->end.dcosx;
  trackInfo->dcosYEnd = ntpTrack->end.dcosy;
  trackInfo->dcosZEnd = ntpTrack->end.dcosz;
  trackInfo->passedFit = 0;
  if(ntpTrack->fit.pass) trackInfo->passedFit = 1;
  trackInfo->reducedChi2 = ntpTrack->fit.chi2;
  if(ntpTrack->fit.ndof>0.) 
    trackInfo->reducedChi2 /= 1.*ntpTrack->fit.ndof;

  trackInfo->forwardRMS = ntpTrack->time.forwardRMS;
  trackInfo->backwardRMS = ntpTrack->time.backwardRMS;

    GetTrackTrace(ntpTrack, trackInfo, -1);
  GetTrackTrace(ntpTrack, trackInfo, 1);

  FillFiducialInformation(trackInfo);

  //loop over the strips in the track and record which planes they are on
  NtpSRStrip *strip = 0;
  for(int i = 0; i < trackInfo->totalStrips; ++i){
    if (ntpTrack->stp[i] >=0) strip = dynamic_cast<NtpSRStrip *>
                                (fStripArray->At(ntpTrack->stp[i]));
    else continue;
    ++trackInfo->stripsPerPlane[strip->plane];
  }
    
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a shower
void ANtpInfoObjectFiller::FillShowerInformation(NtpSRShower *ntpShower,
                                                 ANtpShowerInfo *showerInfo)
{
  showerInfo->Reset();
  showerInfo->planes = ntpShower->plane.n;
  showerInfo->totalStrips = ntpShower->nstrip;
  showerInfo->begPlane = ntpShower->plane.beg;
  showerInfo->endPlane = ntpShower->plane.end;
  showerInfo->vtxX = ntpShower->vtx.x;
  showerInfo->vtxY = ntpShower->vtx.y;
  showerInfo->vtxZ = ntpShower->vtx.z;
  showerInfo->dcosX = ntpShower->vtx.dcosx;
  showerInfo->dcosY = ntpShower->vtx.dcosy;
  showerInfo->dcosZ = ntpShower->vtx.dcosz;
  showerInfo->pulseHeight = ntpShower->ph.sigcor;//used to be sigmap
  showerInfo->linearCCGeV = ntpShower->shwph.linCCgev;
  showerInfo->linearNCGeV = ntpShower->shwph.linNCgev;
  showerInfo->deweightCCGeV = ntpShower->shwph.wtCCgev;
  showerInfo->deweightNCGeV = ntpShower->shwph.wtNCgev;

  //loop over the strips in the shower and record which planes they are on
  NtpSRStrip *strip = 0;
  for(int i = 0; i < showerInfo->totalStrips; ++i){
    strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(ntpShower->stp[i]));
    ++showerInfo->stripsPerPlane[strip->plane];
  }
  
  return;
}

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

  eventInfo->Reset();
  eventInfo->index = static_cast <int> (ntpEvent->index);
  eventInfo->pulseHeight = ntpEvent->ph.sigcor;
  eventInfo->energyGeV = ntpEvent->ph.gev;
  eventInfo->begPlane = ntpEvent->plane.beg;
  eventInfo->endPlane = ntpEvent->plane.end;
  eventInfo->planes = ntpEvent->plane.n;
  
  //add up the pulse height in each plane to find the maximum in the event
  eventInfo->totalStrips = ntpEvent->nstrip;
  
  //write your own criterea for whether a strip is good or not
  eventInfo->passStrips = ntpEvent->nstrip;
  
  eventInfo->showers = ntpEvent->nshower;
  eventInfo->tracks = ntpEvent->ntrack;

  //Correct cedar vertex bug 
  TString reco = ntpManip->GetReleaseMCType();
  if(reco.Contains((ReleaseType::AsString(ReleaseType::kCedarData)).c_str()) || 
     reco.Contains((ReleaseType::AsString(ReleaseType::kCedarCarrot)).c_str()) ||
     reco.Contains((ReleaseType::AsString(ReleaseType::kCedarDaikon)).c_str())){
    
    NtpVtxFinder::NtpVtxFinder vtxf;
    vtxf.SetTargetEvent(static_cast <int> (ntpEvent->index), 
                   ntpManip->GetNtpStRecord());
    if (vtxf.FindVertex()<0){ //Take existing vertex if fix cannot be used
      eventInfo->vtxX = ntpEvent->vtx.x;
      eventInfo->vtxY = ntpEvent->vtx.y;
      eventInfo->vtxZ = ntpEvent->vtx.z;
      eventInfo->vertexTime = ntpEvent->vtx.t;
    } else {
      eventInfo->vtxX = vtxf.VtxX();
      eventInfo->vtxY = vtxf.VtxY();
      eventInfo->vtxZ = vtxf.VtxZ();
      eventInfo->vertexTime = vtxf.VtxT();
    }//if(vtxf.FindVertex())
    
  } else {
    eventInfo->vtxX = ntpEvent->vtx.x;
    eventInfo->vtxY = ntpEvent->vtx.y;
    eventInfo->vtxZ = ntpEvent->vtx.z;
    eventInfo->vertexTime = ntpEvent->vtx.t;
  }//if(reco.Contains(...))

  FillFiducialInformation(eventInfo);

  if(!fStripArray) return;

  //loop over the strips in the shower and record which planes they are on
  NtpSRStrip *strip = 0;
  for(int i = 0; i < eventInfo->totalStrips; ++i){
    if (ntpEvent->stp[i] >= 0) strip = dynamic_cast<NtpSRStrip *>
                                  (fStripArray->At(ntpEvent->stp[i]));
    else continue;
    ++eventInfo->stripsPerPlane[strip->plane];
  }
  
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to the MC truth
void ANtpInfoObjectFiller::FillMCTruthInformation(NtpMCTruth *ntpMCTruth,
                                                  ANtpTruthInfo *truthInfo)
{
  const Float_t muMass = 0.105658357; //mu mass in GeV/c
  const Float_t eMass = 0.000510999; //e mass in GeV/c
  const Float_t tauMass = 1.777;      //tau mass in GeV/c

  truthInfo->Reset();
  truthInfo->nuEnergy = ntpMCTruth->p4neu[3];
  truthInfo->nuVtxX = ntpMCTruth->vtxx;
  truthInfo->nuVtxY = ntpMCTruth->vtxy;
  truthInfo->nuVtxZ = ntpMCTruth->vtxz;
  if(TMath::Abs(truthInfo->nuEnergy) > 0.){
    truthInfo->nuDCosX = ntpMCTruth->p4neu[0]/TMath::Sqrt(truthInfo->nuEnergy*truthInfo->nuEnergy);
    truthInfo->nuDCosY = ntpMCTruth->p4neu[1]/TMath::Sqrt(truthInfo->nuEnergy*truthInfo->nuEnergy);
    truthInfo->nuDCosZ = ntpMCTruth->p4neu[2]/TMath::Sqrt(truthInfo->nuEnergy*truthInfo->nuEnergy);             
  }
  else{
    truthInfo->nuDCosX = -10000.;
    truthInfo->nuDCosY = -10000.;
    truthInfo->nuDCosZ = -10000.;
  }
  truthInfo->nuFlavor = ntpMCTruth->inu;
  truthInfo->interactionType = ntpMCTruth->iaction;
  
  truthInfo->targetEnergy = ntpMCTruth->p4tgt[3];
  truthInfo->targetPX = ntpMCTruth->p4tgt[0];
  truthInfo->targetPY = ntpMCTruth->p4tgt[1];
  truthInfo->targetPZ = ntpMCTruth->p4tgt[2];
  
  truthInfo->hadronicY = ntpMCTruth->y;
  
  truthInfo->showerEnergy = 0.;
  float showerMomentumSqr = 0.;
  if(ntpMCTruth->p4shw[3]>0. && TMath::Abs(ntpMCTruth->p4shw[0] + 
                                           ntpMCTruth->p4shw[1] + 
                                           ntpMCTruth->p4shw[2]) > 0.){
    truthInfo->showerEnergy = ntpMCTruth->p4shw[3];
    showerMomentumSqr = ntpMCTruth->p4shw[0]*ntpMCTruth->p4shw[0];
    showerMomentumSqr += ntpMCTruth->p4shw[1]*ntpMCTruth->p4shw[1];
    showerMomentumSqr += ntpMCTruth->p4shw[2]*ntpMCTruth->p4shw[2];
    truthInfo->showerDCosX = ntpMCTruth->p4shw[0]/TMath::Sqrt(showerMomentumSqr);
    truthInfo->showerDCosY = ntpMCTruth->p4shw[1]/TMath::Sqrt(showerMomentumSqr);
    truthInfo->showerDCosZ = ntpMCTruth->p4shw[2]/TMath::Sqrt(showerMomentumSqr);
  }
  else{
    truthInfo->showerDCosX = ANtpDefVal::kFloat;
    truthInfo->showerDCosY = ANtpDefVal::kFloat;
    truthInfo->showerDCosZ = ANtpDefVal::kFloat;
  }
  
  truthInfo->leptonMomentum = 0.;
  if(TMath::Abs(ntpMCTruth->p4mu1[3])>muMass){
    truthInfo->leptonMomentum = TMath::Sqrt(ntpMCTruth->p4mu1[3]*ntpMCTruth->p4mu1[3] - muMass*muMass);
    truthInfo->leptonDCosX = ntpMCTruth->p4mu1[0]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosY = ntpMCTruth->p4mu1[1]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosZ = ntpMCTruth->p4mu1[2]/truthInfo->leptonMomentum;            
    if(ntpMCTruth->p4mu1[3]<0.) truthInfo->leptonMomentum *= -1.;
  }
  else if(TMath::Abs(ntpMCTruth->p4el1[3])>eMass){
    truthInfo->leptonMomentum = TMath::Sqrt(ntpMCTruth->p4el1[3]*ntpMCTruth->p4el1[3] - eMass*eMass);
    truthInfo->leptonDCosX = ntpMCTruth->p4el1[0]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosY = ntpMCTruth->p4el1[1]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosZ = ntpMCTruth->p4el1[2]/truthInfo->leptonMomentum;            
    if(ntpMCTruth->p4el1[3]<0.) truthInfo->leptonMomentum *= -1.;
  }
  else if(TMath::Abs(ntpMCTruth->p4tau[3])>tauMass){
    truthInfo->leptonMomentum = TMath::Sqrt(ntpMCTruth->p4tau[3]*ntpMCTruth->p4tau[3] - tauMass*tauMass);
    truthInfo->leptonDCosX = ntpMCTruth->p4tau[0]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosY = ntpMCTruth->p4tau[1]/truthInfo->leptonMomentum;
    truthInfo->leptonDCosZ = ntpMCTruth->p4tau[2]/truthInfo->leptonMomentum;            
    if(ntpMCTruth->p4tau[3]<0.) truthInfo->leptonMomentum *= -1.;
  }
  else{
    truthInfo->leptonDCosX = ANtpDefVal::kFloat;
    truthInfo->leptonDCosY = ANtpDefVal::kFloat;
    truthInfo->leptonDCosZ = ANtpDefVal::kFloat;                
  }
    
  
  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFiller::FillFiducialInformation(ANtpTrackInfo *trackInfo)
{
  double x = trackInfo->vtxX; 
  double y = trackInfo->vtxY;
  double z = trackInfo->vtxZ;

  trackInfo->vtxMetersToBeam = MetersToBeam(fDetector, x, y, z); 
  trackInfo->vtxMetersToCoil = MetersToCoil(fDetector, x, y); 
  trackInfo->vtxMetersToCloseEdge = MetersToCloseEdge(fDetector, x,y, 0);
  
  x = trackInfo->endX; 
  y = trackInfo->endY;
  z = trackInfo->endZ;

  trackInfo->endMetersToBeam = MetersToBeam(fDetector, x, y, z);
  trackInfo->endMetersToCoil = MetersToCoil(fDetector, x, y);
  trackInfo->endMetersToCloseEdge = MetersToCloseEdge(fDetector, x, y, 1);
  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFiller::FillFiducialInformation(ANtpEventInfo *eventInfo)
{
  double x = eventInfo->vtxX; 
  double y = eventInfo->vtxY;
  double z = eventInfo->vtxZ;
  
  eventInfo->vtxMetersToBeam = MetersToBeam(fDetector, x, y, z); 
  eventInfo->vtxMetersToCoil = MetersToCoil(fDetector, x, y); 
  eventInfo->vtxMetersToCloseEdge = MetersToCloseEdge(fDetector, x, y, 0);

  return;
}

//----------------------------------------------------------------------
void ANtpInfoObjectFiller::SetStripArray(TClonesArray *strips)
{
  fStripArray = strips;
  return;
}

//----------------------------------------------------------------------
double ANtpInfoObjectFiller::Sqr(double x)
{
  return x*x;
}

//----------------------------------------------------------------------
int ANtpInfoObjectFiller::Sqr(int x)
{
  return x*x;
}

//----------------------------------------------------------------------
float ANtpInfoObjectFiller::Sqr(float x)
{
  return x*x;
}

//------------------------------------------------------------------------
void ANtpInfoObjectFiller::GetTrackTrace(NtpSRTrack *ntpTrack, 
                                         ANtpTrackInfo *trackInfo,
                                         Int_t direction)
{
  
  //get the trace for the track - this is very far detector centric for now,
  //adapted from C. Howcroft's code
  
  //direction tells you whether to trace back from vertex (-1) or forward from end (1)
  
  //side is as labeled below
  /*
             0
          --------
       7 /        \ 1
        /   y|     \
     6 |     |      | 2
       |     ----   |
      5 \      x   / 3
         \        /
          --------
              4
  */  
  double dcos[3] = {direction*ntpTrack->vtx.dcosx,
                    direction*ntpTrack->vtx.dcosy,
                    direction*ntpTrack->vtx.dcosz}; 
  double dcosUV[3] = {direction*ntpTrack->vtx.dcosu,
                      direction*ntpTrack->vtx.dcosv,
                      direction*ntpTrack->vtx.dcosz}; 
  double endPointXY[3] = {ntpTrack->vtx.x, ntpTrack->vtx.y, ntpTrack->vtx.z};
  double endPointUV[3] = {ntpTrack->vtx.u, ntpTrack->vtx.v, ntpTrack->vtx.z};
  
  if(direction > 0){
    endPointXY[0] = ntpTrack->end.x;
    endPointXY[1] = ntpTrack->end.y;
    endPointXY[2] = ntpTrack->end.z;
    endPointUV[0] = ntpTrack->end.u;
    endPointUV[1] = ntpTrack->end.v;
    endPointUV[2] = ntpTrack->end.z;
    dcos[0] = ntpTrack->end.dcosx;
    dcos[1] = ntpTrack->end.dcosy;
    dcos[2] = ntpTrack->end.dcosz;
    dcosUV[0] = ntpTrack->end.dcosu;
    dcosUV[1] = ntpTrack->end.dcosv;
    dcosUV[2] = ntpTrack->end.dcosz;
  }
  
  int side = -1;
  double xzproj[2]={-999., -999.};   
  double entry[3] = {0.};
  double minpath = 99e99;
  double tmppath = 0.;
  double maxDetectorY = 4.;
  double sidelength = kFarSideLength;
  
  if(TMath::IsNaN(endPointUV[0]) || TMath::IsNaN(endPointUV[1]) 
     || TMath::IsNaN(dcosUV[0]) || TMath::IsNaN(dcosUV[1]) 
     || TMath::IsNaN(dcos[2])) return;
  
  if(TMath::Abs(dcos[1])>1.e-6){ 
    // Y - SIDES
    //SIDE 0 
    xzproj[0] = endPointXY[0] + (dcos[0]/dcos[1])*(maxDetectorY - endPointXY[1]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcos[1])*(maxDetectorY - endPointXY[1]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(endPointXY[0] - xzproj[0]) 
              + Sqr(endPointXY[1] - maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath < minpath){
          minpath  = tmppath;
          entry[0] = xzproj[0];
          entry[2] = xzproj[1];
          entry[1] = maxDetectorY;
          side     = 0;
        }                       
      }      
    }//end side 0
    //SIDE 4 
    xzproj[0] = endPointXY[0] + (dcos[0]/dcos[1])*(-maxDetectorY - endPointXY[1]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcos[1])*(-maxDetectorY - endPointXY[1]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(endPointXY[0] - xzproj[0]) 
              + Sqr(endPointXY[1] + maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[0] = xzproj[0];
          entry[2] = xzproj[1];
          entry[1] = -maxDetectorY;
          side     = 4;
        }                       
      }      
    }//end side 4
  }//end Y sides
  if(TMath::Abs(dcos[0])>1e-6){  // X - SIDES
    //SIDE 2 
    xzproj[0] = endPointXY[1] + (dcos[1]/dcos[0])*(maxDetectorY - endPointXY[0]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcos[0])*(maxDetectorY - endPointXY[0]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = (endPointXY[1] - xzproj[0])    * (endPointXY[1] - xzproj[0]) 
              + (endPointXY[0] - maxDetectorY) * (endPointXY[0] - maxDetectorY) 
              + (endPointXY[2] - xzproj[1])    * (endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[1] = xzproj[0];
          entry[2] = xzproj[1];
          entry[0] = maxDetectorY;
          side     = 2;
        }               
      }       
    }//end side 2      
    //SIDE 6 
    xzproj[0] = endPointXY[1] + (dcos[1]/dcos[0])*(-maxDetectorY - endPointXY[0]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcos[0])*(-maxDetectorY - endPointXY[0]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(endPointXY[1] - xzproj[0]) 
              + Sqr(endPointXY[0] + maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[1] = xzproj[0];
          entry[2] = xzproj[1];
          entry[0] = -maxDetectorY;
          side     = 6;
        }
      }                       
    }//end if side 6      
  }//end x sides
  if(TMath::Abs(dcosUV[0])>1.e-6){  // U - SIDES
    //SIDE 1 
    xzproj[0] = endPointUV[1] + (dcosUV[1]/dcosUV[0])*(maxDetectorY - endPointUV[0]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcosUV[0])*(maxDetectorY - endPointUV[0]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(ntpTrack->vtx.v - xzproj[0]) 
              + Sqr(ntpTrack->vtx.u - maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[0] = kInverseSqrt2*(maxDetectorY - xzproj[0]) ;
          entry[2] = xzproj[1];
          entry[1] = kInverseSqrt2*(maxDetectorY + xzproj[0]) ;
          side     = 1;
        }                       
      }      
    }//end if side1 
    //SIDE 5 
    xzproj[0] = endPointUV[1] + (dcosUV[1]/dcosUV[0])*(-maxDetectorY - endPointUV[0]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcosUV[0])*(-maxDetectorY - endPointUV[0]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(ntpTrack->vtx.v - xzproj[0]) 
              + Sqr(ntpTrack->vtx.u + maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[0] = kInverseSqrt2*(-maxDetectorY - xzproj[0]) ;
          entry[2] = xzproj[1];
          entry[1] = kInverseSqrt2*(-maxDetectorY + xzproj[0]) ;
          side     = 5;
        }                             
      }
    }//end if side 5 
  }//end u sides
  if(TMath::Abs(dcosUV[1])>1.e-6){  // V - SIDES
    //SIDE 7
    xzproj[0] = endPointUV[0] + (dcosUV[0]/dcosUV[1])*(maxDetectorY - endPointUV[1]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcosUV[1])*(maxDetectorY - endPointUV[1]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath =  Sqr(ntpTrack->vtx.u - xzproj[0]) 
              + Sqr(ntpTrack->vtx.v - maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt ( tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[0] = kInverseSqrt2*(xzproj[0]-maxDetectorY) ;
          entry[2] = xzproj[1];
          entry[1] = kInverseSqrt2*(xzproj[0]+maxDetectorY) ;
          side     = 7;
        }                       
      }      
    }//end if side7
    //SIDE 3
    xzproj[0] = endPointUV[0] + (dcosUV[0]/dcosUV[1])*(-maxDetectorY - endPointUV[1]);
    xzproj[1] = endPointXY[2] + (dcos[2]/dcosUV[1])*(-maxDetectorY - endPointUV[1]);
    if(TMath::Abs(xzproj[0]) < sidelength && ((xzproj[1] -  endPointXY[2])*dcos[2]) > 0.0){
      tmppath = Sqr(ntpTrack->vtx.u - xzproj[0]) 
              + Sqr(ntpTrack->vtx.v + maxDetectorY) 
              + Sqr(endPointXY[2] - xzproj[1]);
      if(TMath::Finite(tmppath)){
        tmppath = TMath::Sqrt(tmppath);
        if(tmppath<minpath){
          minpath  = tmppath;
          entry[0] = kInverseSqrt2*(xzproj[0]+maxDetectorY) ;
          entry[2] = xzproj[1];
          entry[1] = kInverseSqrt2*(xzproj[0]-maxDetectorY) ;
          side     = 3;
        }                         
      }    
    }//end if side 3    
  }//end v sides
  
  if(direction < 0){
    trackInfo->traceVtx = TMath::Sqrt(Sqr(endPointXY[0]-entry[0]) + Sqr(endPointXY[1]-entry[1])
                                      + Sqr(endPointXY[2]-entry[2]));
    trackInfo->traceVtxZ = entry[2] - endPointXY[2];
    
  } 
  else if(direction > 0){
    trackInfo->traceEnd = TMath::Sqrt(Sqr(ntpTrack->end.x-entry[0]) + Sqr(ntpTrack->end.y-entry[1])
                                      + Sqr(ntpTrack->end.z-entry[2]));
    trackInfo->traceEndZ = entry[2] - ntpTrack->end.z;
    
  } 
  
  return;
}

//----------------------------------------------------------------------
Float_t ANtpInfoObjectFiller::MetersToBeam(Detector::Detector_t det,
                                           Float_t x, Float_t y,
                                           Float_t z)

  // Calculate the distance to the beam spot at this Z position
  // From my log book:
  // Beam is at x = 1.4485m
  // Beam dy/dz = -0.0578 (angle = -3.31deg) to the z-axis 
  // Beam passes through y = 0 at back of golden region (plane 90)
  //  \-> Beam passes through y = 0.3087 @ z = 0  
  // Still looking for hard evidence about this.
{
  if (det == Detector::kNear) {
    const Float_t dydz(-0.0578); 
    const Float_t beamCenterX(1.4885 * Munits::m);
    const Float_t beamCenterY0(0.3087 * Munits::m);

    x -= beamCenterX;
    y -= beamCenterY0 + dydz * z; 
    return TMath::Sqrt(x*x + y*y);
  }

  else if (det == Detector::kFar){
    //Beam center is faintly meaningless for FD cut
    return ANtpDefVal::kFloat;
  }

  else return ANtpDefVal::kFloat;
}
//----------------------------------------------------------------------
Float_t ANtpInfoObjectFiller::MetersToCoil(Detector::Detector_t det,
                                           Float_t x, Float_t y)
{
  if (det == Detector::kNear){
   // Take coil to be at (0,0) Looks about right...  (?) 
    return TMath::Sqrt(x*x + y*y);
  }

  else if (det == Detector::kFar){
    return TMath::Sqrt(x*x + y*y);    
  }

  else return ANtpDefVal::kFloat;
}       
                                   
//----------------------------------------------------------------------
Float_t ANtpInfoObjectFiller::MetersToCloseEdge(Detector::Detector_t det,
                                                Float_t x, Float_t y,
                                                Bool_t NDUseFullPlane)
{
  // Calculate the distance to the nearest (outside) plane edge
  // using the PlaneOutline class. 
  // Currently takes distance to nearest edges seperatly for U and V
  // views, then takes the average of the two distances.
  // It would be better to take a single distance to an average plane
  // but this is not yet possible.

  if (det == Detector::kNear || det == Detector::kFar ){
    
    // Pick a plane coverage
    PlaneCoverage::PlaneCoverage_t outline(PlaneCoverage::kComplete);
    if (det == Detector::kNear){
      if (NDUseFullPlane) outline = PlaneCoverage::kNearFull;
      else outline = PlaneCoverage::kNearPartial;
    }

    PlaneOutline po;  
    
    // Calculate distance (to nearest outside edge) in U. 
    Float_t distU(0.);
    Float_t xedge; Float_t yedge; 
    po.DistanceToOuterEdge(x, y, PlaneView::kU, outline,
                           distU, xedge, yedge);

    // Calculate distance in V.
    Float_t distV(0.);
    po.DistanceToOuterEdge(x, y, PlaneView::kV, outline,
                           distV, xedge, yedge); 

    // Now work out if the point is inside or outside each outline
    // Firstly deal with the coil hole. We consider it inside,
    // but the PlaneOutline class calls it outside.
    
    Bool_t insideU(0.); 
    if ( outline != PlaneCoverage::kNearPartial ) {
      // circle of radius 1/sqrt(2) will enclose the coilHole, 
      // but is itself enclosed by the outer edges
      insideU = ( (x*x + y*y) < 0.5 * Munits::m2 );
      //cout << "X: " << x << ";   Y: " 
      //           << y << ";   inside:" << insideU << endl;
    }
    
    Bool_t insideV(insideU);
    if ( !insideU ){ 
      insideU = po.IsInside(x, y, PlaneView::kU, outline);
      insideV = po.IsInside(x, y, PlaneView::kV, outline);
    }
    distU *= ( insideU ? 1. : -1. );
    distV *= ( insideV ? 1. : -1. );

    //    cout << ">U<: " << insideU << "   " << distU << "\n";
    //cout << ">V<: " << insideV << "   " << distV;
    
    return ( distU + distV ) / 2.;
  }  // if ( @ near or far detectors )

  else return ANtpDefVal::kFloat;    
}

---