AlgShowerSR.cxx

Go to the documentation of this file.

// $Id: AlgShowerSR.cxx,v 1.36 2007/02/04 06:01:05 rhatcher Exp $
//
// AlgShowerSR.cxx
//
// AlgShowerSR is a concrete ShowerSR Algorithm class.
//
// Author:  R. Lee 2001.02.21

#include <cassert>
#include <cmath>

#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
#include "CandShowerSR/AlgShowerSR.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "Conventions/PlaneView.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "RecoBase/CandClusterHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/LinearFit.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "RecoBase/PropagationVelocity.h"

#include "Validity/VldContext.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Navigation/XxxItr.h"

ClassImp(AlgShowerSR)

CVSID("$Id: AlgShowerSR.cxx,v 1.36 2007/02/04 06:01:05 rhatcher Exp $");

//______________________________________________________________________
AlgShowerSR::AlgShowerSR()
{
}

//______________________________________________________________________
AlgShowerSR::~AlgShowerSR()
{
}

//______________________________________________________________________
void AlgShowerSR::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
{

  assert(ch.InheritsFrom("CandShowerHandle"));
  CandShowerHandle &csh = dynamic_cast<CandShowerHandle &>(ch);

  assert(cx.GetDataIn());
  assert(cx.GetDataIn()->InheritsFrom("TObjArray"));

  // grab parameters

  //  Double_t MIPperGeV = ac.GetDouble("MIPperGeV");
  Int_t vtxplanespan = ac.GetInt("VtxPlaneSpan");
  Int_t isCosmic = ac.GetInt("IsCosmic");

  const TObjArray *tary =
                         dynamic_cast<const TObjArray*>(cx.GetDataIn());

  Int_t ubegplane=0;
  Int_t vbegplane=0;
  Int_t uendplane=0;
  Int_t vendplane=0;
  CandStripHandle *begstrip=0;
  CandStripHandle *endstrip=0;

  for (Int_t i=0; i<=tary->GetLast(); i++) {
    TObject *tobj = tary->At(i);
    assert(tobj->InheritsFrom("CandClusterHandle"));
    CandClusterHandle *clusterhandle =
                                 dynamic_cast<CandClusterHandle*>(tobj);
    csh.AddCluster(clusterhandle);
    switch (clusterhandle->GetPlaneView()) {
    case PlaneView::kU:
      if (!ubegplane || clusterhandle->GetBegPlane()<ubegplane) {
        ubegplane = clusterhandle->GetBegPlane();
      }
      if(!uendplane || clusterhandle->GetEndPlane()>uendplane) {
        uendplane = clusterhandle->GetEndPlane();
      }
      break;
    case PlaneView::kV:
      if (!vbegplane || clusterhandle->GetBegPlane()<vbegplane) {
        vbegplane = clusterhandle->GetBegPlane();
      }
      if(!vendplane || clusterhandle->GetEndPlane()>vendplane) {
        vendplane = clusterhandle->GetEndPlane();
      }
      break;
    default:
      break;
    }
    if (!begstrip) {
      CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
      begstrip = stripItr();
    }
  }

  // bail if these shower has an empty daughter list
  if(!begstrip){
    MSG("ShowerSR", Msg::kError)
                 << "Found empty daughter list in CandShowerSR" << endl;
    csh.SetDirCosU(0.);
    csh.SetDirCosV(0.);
    csh.SetDirCosZ(1.);
    csh.SetVtxPlane(-1);
    csh.SetVtxZ(-1.);
    csh.SetVtxU(0.);
    csh.SetVtxV(0.);
    csh.SetVtxT(0.);
    csh.SetEndZ(-1.);
    csh.SetEndU(0.);
    csh.SetEndV(0.);
    csh.SetEndT(0.);
    csh.SetEnergy(0);

    return;
  }

  const VldContext * vld = begstrip->GetVldContext();
  UgliGeomHandle ugh(*vld);

  begstrip = 0;

  Double_t uvtx=0.;
  Double_t vvtx=0.;
  Double_t uvtxph=0.;
  Double_t vvtxph=0.;
  Double_t tvtx=0.;

  Double_t uend=0.;
  Double_t vend=0.;
  Double_t uendph=0.;
  Double_t vendph=0.;
  Double_t tend=0.;

  Bool_t firstvtx(1);
  Bool_t firstend(1);

  Double_t *uvtxfit = new Double_t[vtxplanespan+1];
  Double_t *uvtxphfit = new Double_t[vtxplanespan+1];
  Double_t *uvtxzfit = new Double_t[vtxplanespan+1];
  Double_t *vvtxfit = new Double_t[vtxplanespan+1];
  Double_t *vvtxphfit = new Double_t[vtxplanespan+1];
  Double_t *vvtxzfit = new Double_t[vtxplanespan+1];

  Double_t *uendfit = new Double_t[vtxplanespan+1];
  Double_t *uendphfit = new Double_t[vtxplanespan+1];
  Double_t *uendzfit = new Double_t[vtxplanespan+1];
  Double_t *vendfit = new Double_t[vtxplanespan+1];
  Double_t *vendphfit = new Double_t[vtxplanespan+1];
  Double_t *vendzfit = new Double_t[vtxplanespan+1];

  for (Int_t i=0; i<=vtxplanespan; i++) {
    uvtxfit[i] = 0.;
    uvtxphfit[i] = 0.;
    uvtxzfit[i] = 0.;
    vvtxfit[i] = 0.;
    vvtxphfit[i] = 0.;
    vvtxzfit[i] = 0.;
    uendfit[i] = 0.;
    uendphfit[i] = 0.;
    uendzfit[i] = 0.;
    vendfit[i] = 0.;
    vendphfit[i] = 0.;
    vendzfit[i] = 0.;
  }

  // loop over cluster daughters, calculating a charge sum, and
  // constructing arrays of charge-weighted u vs z and v vs z for the 
  // planes within VtxPlaneSpan of beginning and end of shower

  Double_t charge=0.;
  for (Int_t i=0; i<=tary->GetLast(); i++) {
    TObject *tobj = tary->At(i);
    assert(tobj->InheritsFrom("CandClusterHandle"));
    CandClusterHandle *clusterhandle =
                                 dynamic_cast<CandClusterHandle*>(tobj);
    CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
    while (CandStripHandle *strip = stripItr()) {
      csh.AddDaughterLink(*strip);
      charge += strip->GetCharge();
      if (!begstrip || strip->GetPlane()<begstrip->GetPlane()) {
        begstrip = strip;
      }
      if (!endstrip || strip->GetPlane()>endstrip->GetPlane()) {
        endstrip = strip;
      }

      switch (strip->GetPlaneView()) {
      case PlaneView::kU:
        if (strip->GetPlane()<=ubegplane+vtxplanespan) {
          Double_t charge = strip->GetCharge();
          Double_t tpos = strip->GetTPos();
          uvtx += tpos*charge;
          uvtxph += charge;
          if (firstvtx || strip->GetBegTime()<tvtx) {
            firstvtx = 0;
            tvtx = strip->GetBegTime();
          }
          Int_t dplane = strip->GetPlane()-ubegplane;
          if(dplane<=vtxplanespan && dplane>+0){
            uvtxfit[dplane] += tpos*charge;
            uvtxphfit[dplane] += charge;
            UgliScintPlnHandle upid = 
              ugh.GetScintPlnHandle(strip->GetStripEndId());
            uvtxzfit[dplane] = upid.GetZ0();
          }
        }
        if(strip->GetPlane()>=uendplane-vtxplanespan) {
          Double_t charge = strip->GetCharge();
          Double_t tpos = strip->GetTPos();
          uend += tpos*charge;
          uendph += charge;
          if (firstend || strip->GetBegTime()<tend) {
            firstend = 0;
            tend = strip->GetBegTime();
          }
          Int_t dplane = uendplane-strip->GetPlane();
          if(dplane<=vtxplanespan && dplane>=0){
            uendfit[dplane] += tpos*charge;
            uendphfit[dplane] += charge;
            UgliScintPlnHandle upid = 
              ugh.GetScintPlnHandle(strip->GetStripEndId());
            uendzfit[dplane] = upid.GetZ0();
          }
        }
        break;
      case PlaneView::kV:
        if (strip->GetPlane()<=vbegplane+vtxplanespan) {
          Double_t charge = strip->GetCharge();
          Double_t tpos = strip->GetTPos();
          vvtx += tpos*charge;
          vvtxph += charge;
          if (firstvtx || strip->GetBegTime()<tvtx) {
            firstvtx = 0;
            tvtx = strip->GetBegTime();
          }
          Int_t dplane = strip->GetPlane()-vbegplane;
          if(dplane<=vtxplanespan && dplane>=0){
            vvtxfit[dplane] += tpos*charge;
            vvtxphfit[dplane] += charge;
            UgliScintPlnHandle upid = 
              ugh.GetScintPlnHandle(strip->GetStripEndId());
            vvtxzfit[dplane] = upid.GetZ0();
          }
        }
        if(strip->GetPlane()>=uendplane-vtxplanespan) {
          Double_t charge = strip->GetCharge();
          Double_t tpos = strip->GetTPos();
          vend += tpos*charge;
          vendph += charge;
          if (firstend || strip->GetBegTime()<tend) {
            firstend = 0;
            tend = strip->GetBegTime();
          }
          Int_t dplane = vendplane-strip->GetPlane();
          if(dplane<=vtxplanespan && dplane>=0){
            vendfit[dplane] += tpos*charge;
            vendphfit[dplane] += charge;
            UgliScintPlnHandle upid = 
              ugh.GetScintPlnHandle(strip->GetStripEndId());
            vendzfit[dplane] = upid.GetZ0();
          }
        }
        break;
      default:
        break;
      }
    }
  }

  // normalize charge weighted fit points, and remove array entries
  // with zero pulse height

  Int_t iuvtxfit=0;
  Int_t ivvtxfit=0;
  Int_t iuendfit=0;
  Int_t ivendfit=0;

  for (Int_t i=0; i<=vtxplanespan; i++) {
    if (uvtxphfit[i]>0.) {
      uvtxfit[iuvtxfit] = uvtxfit[i]/uvtxphfit[i];
      uvtxzfit[iuvtxfit]=uvtxzfit[i];
      uvtxphfit[iuvtxfit]=uvtxphfit[i];

      iuvtxfit++;
    }
    if (vvtxphfit[i]>0.) {
      vvtxfit[ivvtxfit] = vvtxfit[i]/vvtxphfit[i];
      vvtxzfit[ivvtxfit]=vvtxzfit[i];
      vvtxphfit[ivvtxfit]=vvtxphfit[i];
      ivvtxfit++;
    }
   if (uendphfit[i]>0.) {
      uendfit[iuendfit] = uendfit[i]/uendphfit[i];
      uendzfit[iuendfit]=uendzfit[i];
      uendphfit[iuendfit]=uendphfit[i];
      iuendfit++;
    }
    if (vendphfit[i]>0.) {
      vendfit[ivendfit] = vendfit[i]/vendphfit[i];
      vendzfit[ivendfit]=vendzfit[i];
      vendphfit[ivendfit]=vendphfit[i];
      ivendfit++;
    }
  }

  //perform weighted fit of u, v vs z positions
  Double_t uvtxparm[2];
  Double_t vvtxparm[2];
  Double_t dudz=0.;
  Double_t dvdz=0.;
  if (iuvtxfit>1) {
    LinearFit::Weighted(iuvtxfit,uvtxzfit,uvtxfit,uvtxphfit,uvtxparm);
    dudz = uvtxparm[1];
  }
  if (ivvtxfit>1) {
    LinearFit::Weighted(ivvtxfit,vvtxzfit,vvtxfit,vvtxphfit,vvtxparm);
    dvdz = vvtxparm[1];
  }
  Double_t uendparm[2];
  Double_t vendparm[2];
  if (iuendfit>1) {
    LinearFit::Weighted(iuendfit,uendzfit,uendfit,uendphfit,uendparm);
 
  }
  if (ivendfit>1) {
    LinearFit::Weighted(ivendfit,vendzfit,vendfit,vendphfit,vendparm);
 
  }

  // calculate direction cosines from linear fits above
  Double_t cosz=1./sqrt(1.+dudz*dudz+dvdz*dvdz);
  Double_t cosu=dudz*cosz;
  Double_t cosv=dvdz*cosz;

  csh.SetDirCosU(cosu);
  csh.SetDirCosV(cosv);
  csh.SetDirCosZ(cosz);

  csh.SetEndDirCosU(cosu);
  csh.SetEndDirCosV(cosv);
  csh.SetEndDirCosZ(cosz);

 // set the vertex planes and z position.  The vertex Z position
  // is the position of the most upstream strip.  This is not
  // correct for cosmics!

  PlexStripEndId stripid = begstrip->GetStripEndId();
  UgliScintPlnHandle planehandle = ugh.GetScintPlnHandle(stripid);
  csh.SetVtxPlane(begstrip->GetPlane());
  csh.SetVtxZ(planehandle.GetZ0());

  stripid = endstrip->GetStripEndId();
  planehandle = ugh.GetScintPlnHandle(stripid);
  csh.SetEndPlane(endstrip->GetPlane());
  csh.SetEndZ(planehandle.GetZ0());

  // the vertex time is defined to be the earliest strip begin time in
  // the vertex region defined by VtxPlaneSpan
  csh.SetVtxT(tvtx);
  csh.SetEndT(tend);  
  SetUV(&csh);
 
  // if this is a cosmic event, base direction on timing
  if(isCosmic){
    Double_t fitparm[2] = {0,0};
    Double_t timefitchi2;
    FindTimingDirection(csh,
                        ac,
                        fitparm,
                        timefitchi2);
    if(fitparm[1]<0){
      
      MSG("AlgShowerSR",Msg::kDebug) << " swapping vertex " << endl;

      Int_t newvtx = csh.GetEndPlane();
      Int_t newend = csh.GetVtxPlane();
      Double_t newvtxz = csh.GetEndZ();
      Double_t newendz = csh.GetVtxZ();
      Double_t newvtxt = csh.GetEndT();
      Double_t newendt = csh.GetVtxT();

      csh.SetVtxZ(newvtxz);
      csh.SetEndZ(newendz);
      csh.SetVtxPlane(newvtx);
      csh.SetEndPlane(newend);
      csh.SetVtxT(newvtxt);
      csh.SetEndT(newendt);  
    }
  }

  
  // the vertex and end U/V positions are based on fits if the fit was done

  csh.SetVtxU(csh.GetU(csh.GetVtxPlane()));
  csh.SetVtxV(csh.GetV(csh.GetVtxPlane()));
  
  csh.SetEndU(csh.GetU(csh.GetEndPlane()));
  csh.SetEndV(csh.GetV(csh.GetEndPlane()));

  Calibrate(&csh);
  
  CandTrackHandle * associatedtrk=0;
  csh.CalibrateEnergy(associatedtrk,ac);

  delete [] uvtxfit;
  delete [] uvtxphfit;
  delete [] uvtxzfit;
  delete [] vvtxfit;
  delete [] vvtxphfit;
  delete [] vvtxzfit;

  delete [] uendfit;
  delete [] uendphfit;
  delete [] uendzfit;
  delete [] vendfit;
  delete [] vendphfit;
  delete [] vendzfit;
}

//______________________________________________________________________

void AlgShowerSR::SetUVT(CandShowerHandle *candshower)
{
  this->SetUV(candshower);
  this->SetT(candshower);
}

//______________________________________________________________________

void AlgShowerSR::SetUV(CandShowerHandle *candshower)
{
  TIter stripItr(candshower->GetDaughterIterator());
  CandStripHandle *firststrip = 0;
  firststrip = dynamic_cast<CandStripHandle*>(stripItr());
  if (!firststrip) return;  // if no strips, do nothing
  
  const VldContext *vldc = firststrip->GetVldContext();
  Detector::Detector_t det = vldc->GetDetector();

  UgliGeomHandle ugh = UgliGeomHandle(*vldc);

  CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator());
  CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
  stripKf->SetFun(CandStripHandle::KeyFromPlane);
  orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
  stripKf = 0;
  Bool_t first(1);
  Int_t oldplane = -1;
  PlaneView::PlaneView_t oldview = PlaneView::kUnknown;
  Float_t tpos = 0.;
  Float_t ph = 0.;
  Int_t idir=1;
  Int_t vtxplane = candshower->GetBegPlane();
  Int_t endplane = candshower->GetEndPlane();
  if (vtxplane>endplane) idir=-1;
  Int_t begplane = min(vtxplane,endplane);
  endplane = max(vtxplane,endplane);

  begplane = -1;
  endplane = -1;

// find begin and end planes
  while (CandStripHandle *strip =
                    dynamic_cast<CandStripHandle*>(orderedstripItr())) {
    if (strip->GetPlaneView()==PlaneView::kU ||
        strip->GetPlaneView()==PlaneView::kV) {
      if (begplane<0 || strip->GetPlane()<begplane) {
        begplane = strip->GetPlane();
      }
      if (endplane<0 || strip->GetPlane()>endplane) {
        endplane = strip->GetPlane();
      }
    }
  }
  Int_t *planelist = new Int_t[endplane+1];
  Int_t *stripendlist = new Int_t[endplane+1];
  if(begplane>-1 && endplane>-1){

    for (int i=0; i<=endplane; i++) {
      planelist[i] = 0;
      stripendlist[i] = 0;
    }
    
    // transverse positions
    orderedstripItr.Reset();
    
    while (CandStripHandle *strip =
           dynamic_cast<CandStripHandle*>(orderedstripItr())) {
      if (strip->GetPlaneView()==PlaneView::kU ||
          strip->GetPlaneView()==PlaneView::kV) {
        planelist[strip->GetPlane()] = strip->GetPlaneView();
        if (strip->GetNDigit(StripEnd::kNegative)) {
          stripendlist[strip->GetPlane()] += (1<<StripEnd::kNegative);
        }
        if (strip->GetNDigit(StripEnd::kPositive)) {
          stripendlist[strip->GetPlane()] += (1<<StripEnd::kPositive);
        }
        if (first || strip->GetPlane()==oldplane) {
          if (first) {
            first = 0;
            oldplane = strip->GetPlane();
            oldview = strip->GetPlaneView();
          }
          Float_t stripph = strip->GetCharge();
          tpos += strip->GetTPos()*stripph;
          ph += stripph;
        }
        else {
          tpos /= ph;
          if (oldview==PlaneView::kU) {
            candshower->SetU(oldplane,tpos);
          } else if (oldview==PlaneView::kV) {
            candshower->SetV(oldplane,tpos);
          }
          oldplane = strip->GetPlane();
          oldview = strip->GetPlaneView();
          Float_t stripph = strip->GetCharge();
          tpos = strip->GetTPos()*stripph;
          ph = stripph;
        }
      }
    }
    if (ph>0.) {
      tpos /= ph;
      if (oldview==PlaneView::kU) {
        candshower->SetU(oldplane,tpos);
      } else if (oldview==PlaneView::kV) {
        candshower->SetV(oldplane,tpos);
      }
    }
    
    // create a starting and stopping PlexPlaneId
    Int_t idirAdjoin = (begplane>endplane) ? -1 : 1; // direction
    PlexPlaneId planeid(det,begplane,false);
    PlexPlaneId planeidend(det,endplane,false);
    PlexPlaneId planeid_toofar = planeidend.GetAdjoinScint(idirAdjoin);
    
    // after each step, if not reached end, move to adjoining *scint* plane
    for ( ; planeid != planeid_toofar ; 
          planeid  = planeid.GetAdjoinScint(idirAdjoin) ) {
      
      // quit if hit something meaningless
      if ( ! planeid.IsValid() ) break;
      
      UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(planeid);
      
      if ( ! scintpln.IsValid() ) continue;
      
      Int_t iplane = planeid.GetPlane();
      
      Float_t zpos = scintpln.GetZ0();
      for (Int_t planeview=PlaneView::kU; planeview<=PlaneView::kV;
           planeview++) {
        if (planelist[iplane]!=planeview) {
          Int_t found = 0;
          Int_t nfound = 0;
          Int_t jplane[2];
          
          // search upstream
          for (Int_t iplane1=iplane-1; iplane1>=begplane && !found;
               iplane1--) {
            if (planelist[iplane1]==planeview) {
              found++;
              jplane[nfound] = iplane1;
              nfound++;
            }
          }
          Int_t ndown = 1;
          if (!found) {
            ndown++;
          }
          found = 0;
          
          // search downstream
          for (Int_t iplane1=iplane+1; iplane1<=endplane && found<ndown;
               iplane1++) {
            if (planelist[iplane1]==planeview) {
              found++;
              jplane[nfound] = iplane1;
              nfound++;
            }
          }
          if (!found) {
            for (Int_t iplane1=jplane[0]-1; iplane1>=begplane && !found;
                 iplane1--) {
              if (planelist[iplane1]==planeview) {
                found++;
                jplane[nfound] = iplane1;
                nfound++;
              }
            }
          }
          if (nfound==1) {
            if (planeview==PlaneView::kU) {
              candshower->SetU(iplane,candshower->GetU(jplane[0]));
            } else {
              candshower->SetV(iplane,candshower->GetV(jplane[0]));
            }
          }
          else {
            Float_t z[2];
            z[0] = candshower->GetZ(jplane[0]);
            z[1] = candshower->GetZ(jplane[1]);
            Float_t x[2];
            if (planeview==PlaneView::kU) {
              x[0] = candshower->GetU(jplane[0]);
              x[1] = candshower->GetU(jplane[1]);
              Float_t denom = (z[1]-z[0]);
              if (denom) candshower->SetU(iplane,
                                          x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom));
              else {
                MSG("Shower",Msg::kWarning) << endl
                                            << "Denominator (z[1]-z[0]) is zero."
                                            << "  SetU not called." << endl;
              }
            } else {
              x[0] = candshower->GetV(jplane[0]);
              x[1] = candshower->GetV(jplane[1]);
              Float_t denom = (z[1]-z[0]);
              if (denom) candshower->SetV(iplane,
                                          x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom));
              else {
                MSG("Shower",Msg::kWarning) << endl
                                            << "Denominator is (z[1]-z[0]) zero."
                                            << "  SetV not called." << endl;
              }
            }
          }
        }
      }
    }
  }
  delete [] planelist;
  delete [] stripendlist;

}

//______________________________________________________________________

void AlgShowerSR::SetT(CandShowerHandle *candshower)
{

  // we take a weighted average of hits in the same
  // plane.  the proper way to do this would be to consider hits coming
  // from the same PMT and use the earliest time.

  MSG("Shower",Msg::kDebug) << "starting SetT" << endl;

  TIter stripItr(candshower->GetDaughterIterator());
  CandStripHandle *firststrip = 0;
  firststrip = dynamic_cast<CandStripHandle*>(stripItr());
  UgliGeomHandle *ugh = 0;
  SimFlag::SimFlag_t simFlag = SimFlag::kUnknown;
  if (firststrip) {
    ugh = new UgliGeomHandle(*firststrip->GetVldContext());
    simFlag = firststrip->GetVldContext()->GetSimFlag();
  }

  Double_t timewgt[2] = {0.,0.};
  Double_t time[2] = {0.,0.};
  StripEnd::StripEnd_t stripend[2] = {StripEnd::kNegative,StripEnd::kPositive};

  Double_t apos;
  ArrivalTime arrtime;
  CalTimeType::CalTimeType_t caltimetype = CalTimeType::kNone;
  Double_t propagation_velocity = PropagationVelocity::Velocity(simFlag)*Munits::ns;

  CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator());
  CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
  stripKf->SetFun(CandStripHandle::KeyFromPlane);
  orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
  stripKf = 0;

  Int_t oldplane = -1;
  PlaneView::PlaneView_t oldview = PlaneView::kUnknown;

  Bool_t first(1);
  while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) {
    if (strip->GetPlaneView()==PlaneView::kU ||
        strip->GetPlaneView()==PlaneView::kV) {
      if (strip->GetPlaneView()==PlaneView::kU) {
        apos = candshower->GetV(strip->GetPlane());
      } else {
        apos = candshower->GetU(strip->GetPlane());
      }
      apos=min(apos,4.0); 
      apos=max(apos,-4.0);  // sanity check on apos bounds
      UgliStripHandle striphandle = ugh->GetStripHandle(strip->GetStripEndId());
      caltimetype = strip->GetCalTimeType();
      if (first || strip->GetPlane()==oldplane) {
        if (first) {
          first = 0;
          oldplane = strip->GetPlane();
          oldview = strip->GetPlaneView();
        }
        Double_t thisph[2] = {
          strip->GetCharge(CalDigitType::kPE,stripend[0]),
          strip->GetCharge(CalDigitType::kPE,stripend[1])};
        Double_t thiswgt[2] = {
          min(0.4,arrtime.Weight(max(1.,thisph[0]))),
          min(0.4,arrtime.Weight(max(1.,thisph[1])))};
        for (int i=0; i<2; i++) {
          if (thisph[i]>0.) {
            Double_t plusminus = 1.;
            if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
                (strip->GetPlaneView()==PlaneView::kU && i==0)) {
              plusminus = -1.;
            }
            Double_t thistime = strip->GetBegTime(stripend[i]) - 
              1.e-9*(striphandle.ClearFiber(stripend[i]) + 
                     striphandle.WlsPigtail(stripend[i]) + 
                     (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity;

            MSG("Shower",Msg::kDebug) 
              << "plane " << strip->GetPlane() << " end " << i 
              << " raw time " << strip->GetBegTime(stripend[i])*1e9 
              << " pulse height " << thisph[i] << " weight " << thiswgt[i] 
              << " clearfiber " << striphandle.ClearFiber(stripend[i]) 
              << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) 
              << " halflength " << striphandle.GetHalfLength() 
              << " apos " << apos << " 3D time " << thistime*1e9 << endl;
            time[i] += thiswgt[i]*thistime;
            timewgt[i] += thiswgt[i];
            MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
          }
        }
      }
      else {
        for (int i=0; i<2; i++) {
          if (timewgt[i]>0.) {
            time[i] /= timewgt[i];
            candshower->SetT(oldplane,stripend[i],time[i]);
            MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " 
                                      << stripend[i] << " " << time[i]*1e9 << endl;
          }
        }
        oldplane = strip->GetPlane();
        oldview = strip->GetPlaneView();
        Double_t thisph[2] = {
          strip->GetCharge(CalDigitType::kPE,stripend[0]),
          strip->GetCharge(CalDigitType::kPE,stripend[1])};
        Double_t thiswgt[2] = {
          min(0.4,arrtime.Weight(max(1.,thisph[0]))),
          min(0.4,arrtime.Weight(max(1.,thisph[1])))};
        for (int i=0; i<2; i++) {
          time[i] = 0.;
          timewgt[i] = 0.;
          if (thisph[i]>0.) {
            Double_t plusminus = 1.;
            if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
                (strip->GetPlaneView()==PlaneView::kU && i==0)) {
              plusminus = -1.;
            }
            Double_t thistime = strip->GetBegTime(stripend[i]) - 
              1.e-9*(striphandle.ClearFiber(stripend[i]) + 
                     striphandle.WlsPigtail(stripend[i]) + 
                     (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity;
            MSG("Shower",Msg::kDebug)
              << "plane " << strip->GetPlane() << " end " << i 
              << " raw time " << strip->GetBegTime(stripend[i])*1e9 
              << " pulse height " << thisph[i] << " weight " << thiswgt[i] 
              << " clearfiber " << striphandle.ClearFiber(stripend[i]) 
              << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) 
              << " halflength " << striphandle.GetHalfLength() 
              << " apos " << apos << " 3D time " << thistime*1e9 << endl;
            time[i] = thiswgt[i]*thistime;
            timewgt[i] = thiswgt[i];
            MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
          }
        }
      }
    }
  }
  for (int i=0; i<2; i++) {
    if (timewgt[i]>0.) {
      time[i] /= timewgt[i];
      candshower->SetT(oldplane,stripend[i],time[i]);
      MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " 
                                << stripend[i] << " " << time[i]*1e9 << endl;
    }
  }
  if (ugh) delete ugh;
}

//---------------------------------------------------------------------
//figure out if the shower is going forward or backwards with respect to
//z.  return the number of points used to determine the direction
Int_t AlgShowerSR::FindTimingDirection(CandShowerHandle &shwh,
                                       AlgConfig &ac,
                                       Double_t *fitparm,
                                       Double_t &timefitchi2)
{
   Double_t parmMaxTimeFitRes = ac.GetDouble("MaxTimeFitRes");
  Double_t parmMinTimeFitSize = ac.GetInt("MinTimeFitSize");
  Double_t parmMaxTimeFitOutlier = ac.GetDouble("MaxTimeFitOutlier");
  
  MSG("AlgShowerSR", Msg::kVerbose) << "in find timing direction" << endl;
  
  fitparm[0] = 0;
  fitparm[1] = 0;
  
  CandStripHandle *strip = 0;
  CandDigitHandle *digit = 0;
  
  StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown;
  
  //get VldContext and UgliGeomHandle
  CandStripHandleItr stripItr(shwh.GetDaughterIterator());
  CandStripHandle *firststrip = stripItr.Ptr();
  
  assert(firststrip);
  
  const VldContext *vldcontext = firststrip->GetVldContext();
  assert(vldcontext);
  
  UgliGeomHandle ugh(*vldcontext);
  
  //arrays for determining the timing direction - limit to
  //1000 entries - really shouldnt need more than that to 
  //be able to figure out if the shower is going north or 
  //south
  Double_t time[1000];
  Double_t pathLength[1000];
  Double_t weight[1000];
  Double_t adcsum[1000];
  
  for(Int_t i = 0; i<1000; ++i){
    time[i] = 0.;
    pathLength[i] = 0.;
    weight[i] = 0.;
    adcsum[i] = 0;
  }
  
  //need to put in a check to make sure you dont have more planes in one
  //view than the other
  Int_t uBeg = 1000;
  Int_t uEnd = 0;
  Int_t vBeg = 1000;
  Int_t vEnd = 0;
   while( (strip = stripItr())){
    Int_t plane = strip->GetPlane();
    if(strip->GetPlaneView() == PlaneView::kV && plane>vEnd) vEnd=plane;
    if(strip->GetPlaneView() == PlaneView::kV && plane<vBeg) vBeg=plane;
    if(strip->GetPlaneView() == PlaneView::kU && plane>uEnd) uEnd=plane;
    if(strip->GetPlaneView() == PlaneView::kU && plane<uBeg) uBeg=plane;
  }
  
  //get the endpoints  
  Int_t showerBeg = TMath::Min(uBeg,vBeg);
  Int_t showerEnd = TMath::Max(uEnd,vEnd);
  
  if(TMath::Abs(uBeg-vBeg)>3.){
    showerBeg = TMath::Max(uBeg,vBeg) - 1;
  } 
  if(TMath::Abs(uEnd-vEnd)>3.){
    showerEnd = TMath::Min(uEnd,vEnd) + 1;
  }
  MSG("AlgShowerSR", Msg::kDebug) << "u end points = " << uBeg
                                  << " " << uEnd << " v end points = "
                                  << vBeg << " " << vEnd << endl
                                  << " shower end points = " << showerBeg
                                  << " " << showerEnd << endl;
  Double_t halfLength = 0.;
  Double_t apos = 0.;
  Double_t distFromCenter = 0.;
  Double_t fiberDist = 0.;
  Int_t nctr = 0;
  Double_t maxPathLength = 0;
  for(Int_t iplane = showerBeg;iplane <= showerEnd;iplane++){
    adcsum[nctr] = 0;
    time[nctr] = 0; 
    pathLength[nctr] = 0;
    CandStripHandleItr shwstripItr(shwh.GetDaughterIterator());
    while( (strip = shwstripItr()) && nctr<1000){
      if(strip->GetPlane() == iplane){
        MSG("AlgShowerSR", Msg::kDebug) << "strip from plane " << iplane
                                        << " " << Detector::AsString(vldcontext->GetDetector()) << endl;
        
        //only look at double ended strips if in the far detector
        //and strips from planes that have orthogonal measurements 
        if( (strip->GetNDaughters() == 2 || vldcontext->GetDetector() == Detector::kNear)){
          
          MSG("AlgShowerSR", Msg::kVerbose) << "strip good for timing" << endl;
          
          //get the iterator over the digits in the strip
          CandDigitHandleItr digitItr(strip->GetDaughterIterator());
          
          while( (digit = digitItr())){     
            stripEnd = digit->GetPlexSEIdAltL().GetEnd();
            pathLength[nctr] = strip->GetZPos();
            maxPathLength = TMath::Max(maxPathLength,pathLength[nctr] );
            UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId());
            halfLength = stripHandle.GetHalfLength();
            apos = 0.;
            fiberDist = 0.;
            
            if(strip->GetPlaneView() == PlaneView::kV) apos = shwh.GetU(iplane);
            else if(strip->GetPlaneView() == PlaneView::kU) apos = shwh.GetV(iplane);
            distFromCenter = 0.;
            if(strip->GetPlaneView() == PlaneView::kV) 
              distFromCenter = (stripEnd == StripEnd::kNegative) ? apos : -apos;
            else if(strip->GetPlaneView() == PlaneView::kU) 
              distFromCenter = (stripEnd == StripEnd::kNegative) ? -apos : apos;
            
            fiberDist = (halfLength + distFromCenter + stripHandle.ClearFiber(stripEnd) 
                         + stripHandle.WlsPigtail(stripEnd));
            
            time[nctr] += digit->GetCharge() * (strip->GetTime(stripEnd) - fiberDist/PropagationVelocity::Velocity());   
            
            weight[nctr] = 1.0;
            adcsum[nctr] += digit->GetCharge();
            Double_t temptime =  time[nctr];
            if(adcsum[nctr]!=0)temptime /= adcsum[nctr]; 
            MSG("AlgShowerSR", Msg::kDebug) << nctr 
                                            << " " << pathLength[nctr]
                                            << " " << temptime
                                            << " " << strip->GetTime(stripEnd)
                                            << " " << weight[nctr]
                                            << endl;
            
          }//end loop over digits in strip
        }//end if double ended strip
      }//end loop over strips in plane
    }
    if(adcsum[nctr] > 0){
      time[nctr] /= adcsum[nctr];        
      nctr++;
    } 
  }
  for(Int_t i=0;i<nctr;i++){
    MSG("AlgShowerSR",Msg::kDebug) << "TIMEFIT " << i << " " 
                                   <<pathLength[i] << " " 
                                   << (time[i]-time[0])*1.e9 << " " << weight[i] << endl;
  }
  
  //loop over all strips in the shower
  Double_t efitparm[2];
  Double_t maxresid = parmMaxTimeFitRes + 1.;
  Double_t resid = 0.;
  Int_t ctr = nctr;
  Int_t imaxresid = 0;
  Int_t nremove = 0;
  Bool_t dofit = false;
  while(maxresid > parmMaxTimeFitRes
        && nctr-nremove > parmMinTimeFitSize 
        && (1.*nremove) < parmMaxTimeFitOutlier*nctr
        ){
    dofit = true;
    LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm);
    maxresid = 0.;
    imaxresid = 0;
    for(int i=0; i<ctr; ++i){
      resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9;
      if(weight[i]>0. && TMath::Abs(resid)>maxresid){
        maxresid = TMath::Abs(resid);
        imaxresid = i;
      }
    }
    
    MSG("AlgShowerSR",Msg::kDebug) << "constrained fit, dt/ds slope, offset = " 
                                   << fitparm[1] << " " << fitparm[0] << endl;
    MSG("AlgShowerSR",Msg::kDebug) << "maximum residual " << maxresid << " " 
                                   << pathLength[imaxresid] << " " 
                                   << time[imaxresid] << " " 
                                   << weight[imaxresid] << endl;
    
    if(maxresid>parmMaxTimeFitRes){
      MSG("AlgShowerSR",Msg::kDebug) << "  removing" << endl;
      weight[imaxresid] = 0.;
      nremove++;
    }
    
    timefitchi2 = 0.;
    for(int i = 0; i < ctr; ++i){
      resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1e9;
      timefitchi2 += weight[i]*resid*resid;
    }
    MSG("AlgShowerSR",Msg::kDebug) << "chi2 = " << timefitchi2 \
                                   << "   n = " << nctr-nremove << endl;
  }//end loop to find timing fit and remove outliers
  
  timefitchi2 = 0.;
  if(dofit){
    for(int i = 0; i<ctr; ++i){
      if(weight[i] > 0.){
        resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9;
        timefitchi2 += weight[i]*resid*resid;
      }
    }
  }

  MSG("AlgShowerSR",Msg::kDebug) << " TIMEFIT " << fitparm[1] 
                                 << " chi^2/ndf = " 
                                 << timefitchi2/(1.*(nctr-nremove)) 
                                 << " " << TMath::Abs(uBeg-uEnd)
                                 << " " << TMath::Abs(vBeg-vEnd) << endl;
  
  //check the chi^2 for the fit - if it is really bad dont change the
  //initial guess at the direction for the event, ie just make sure
  //that the slope in time vs pathlength is positive;
  if(fitparm[1] < 0. && timefitchi2 >= 10.*(nctr-nremove))
    fitparm[1] *= -1.;
  return nctr-nremove;
}


//______________________________________________________________________
void AlgShowerSR::Trace(const char * /* c */) const
{
}

---