CandShowerSRHandle.cxx

Go to the documentation of this file.

// $Id: CandShowerSRHandle.cxx,v 1.26 2007/02/04 06:01:05 rhatcher Exp $
//
// CandShowerSRHandle.cxx
//
// CandShowerSRHandle is the specialized access handle to CandShowerSR.
//

#include <cassert>
#include <iostream>

#include "MessageService/MsgService.h"
#include "CandShowerSR/CandShowerSR.h"
#include "CandShowerSR/CandShowerSRHandle.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Validity/VldContext.h"
#include "Algorithm/AlgConfig.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandFitTrackHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "Conventions/Mphysical.h"
#include "Conventions/Munits.h"
#include "TMath.h"

ClassImp(CandShowerSRHandle)

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

//______________________________________________________________________
CandShowerSRHandle::CandShowerSRHandle()
{
}

//______________________________________________________________________
CandShowerSRHandle::CandShowerSRHandle(const CandShowerSRHandle &cdh) :
  CandShowerHandle(cdh)
{
}

//______________________________________________________________________
CandShowerSRHandle::CandShowerSRHandle(CandShowerSR *cd) :
  CandShowerHandle(cd)
{
}

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

//______________________________________________________________________
CandShowerSRHandle *CandShowerSRHandle::DupHandle() const
{
   return (new CandShowerSRHandle(*this));
}

//______________________________________________________________________
void CandShowerSRHandle::Trace(const char *c) const
{
  MSG("Cand", Msg::kDebug)
    << "**********Begin CandShowerSRHandle::Trace(\"" << c << "\")" << endl
           << "Information from CandShowerSRHandle's CandHandle: " << endl;
  CandHandle::Trace(c);
  MSG("Cand", Msg::kDebug)
     << "**********End CandShowerSRHandle::Trace(\"" << c << "\")" << endl;
}

//_____________________________________________________________________
void CandShowerSRHandle::AddSubShower(CandSubShowerSRHandle *subshower)
{
  const TObjArray &subshowerlist =
    (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList;
  Bool_t found(0);
  for (Int_t i=0; i<=subshowerlist.GetLast() && !found; i++) {
    CandSubShowerSRHandle *target =
      dynamic_cast<CandSubShowerSRHandle*>(subshowerlist.At(i));
    if (*subshower == *target) found = 1;
  }
  if (!found) {                 // don't want to duplicate object in list
    if(subshower->GetPlaneView()==PlaneView::kU ||
       subshower->GetPlaneView()==PlaneView::kX) {
      (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList.
        Add(subshower->DupHandle());
      dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nUSubShowers+=1;
    }
    else if(subshower->GetPlaneView()==PlaneView::kV ||
            subshower->GetPlaneView()==PlaneView::kY) {
      (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList.
        Add(subshower->DupHandle());
      dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nVSubShowers+=1;
    }
    else {
      MSG("SubShowerSR",Msg::kWarning) << "Unknown SubShower View! Not adding it to ClusterList"
                                       << endl;     
    }
  }
  return;
}

//______________________________________________________________________
void CandShowerSRHandle::RemoveSubShower(CandSubShowerSRHandle *subshower)
{
  const TObjArray &subshowerlist =
    (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList;
  Bool_t found(0);
  for (Int_t i=0; i<=subshowerlist.GetLast() && !found; i++) {
    CandSubShowerSRHandle *target =
      dynamic_cast<CandSubShowerSRHandle*>(subshowerlist.At(i));
    if (*subshower == *target){
      if(subshower->GetPlaneView()==PlaneView::kU ||
         subshower->GetPlaneView()==PlaneView::kX) {
        dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nUSubShowers-=1;
      }
      else if(subshower->GetPlaneView()==PlaneView::kV ||
              subshower->GetPlaneView()==PlaneView::kY) {
        dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nVSubShowers-=1;
      }
      (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList.
        RemoveAt(i);
      delete target;
      (dynamic_cast<CandShowerSR*>(GetOwnedCandBase()))->fSubShowerList.
        Compress();
      return;
    }
  }
  return;
}

//______________________________________________________________________

const CandSubShowerSRHandle *CandShowerSRHandle::GetSubShower(Int_t i) const
{
  const TObjArray *fSubShowerList = &(dynamic_cast<const CandShowerSR*>
    (GetCandBase())->fSubShowerList);
  if (i>fSubShowerList->GetLast()) {
    return 0;
  }
  return dynamic_cast<const CandSubShowerSRHandle*>(fSubShowerList->At(i));
}

Int_t CandShowerSRHandle::GetLastSubShower() const
{
  return dynamic_cast<const CandShowerSR*>(GetCandBase())->fSubShowerList.
                                                              GetLast();
}

//______________________________________________________________________

void CandShowerSRHandle::SetNumSubShowersU(Int_t nusubshw)
{
  dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nUSubShowers = nusubshw;
}

Int_t CandShowerSRHandle::GetNumSubShowersU() const
{
  return dynamic_cast<const CandShowerSR*>(GetCandBase())->nUSubShowers;
}

//______________________________________________________________________

void CandShowerSRHandle::SetNumSubShowersV(Int_t nusubshw)
{
  dynamic_cast<CandShowerSR*>(GetOwnedCandBase())->nVSubShowers = nusubshw;
}

Int_t CandShowerSRHandle::GetNumSubShowersV() const
{
  return dynamic_cast<const CandShowerSR*>(GetCandBase())->nVSubShowers;
}

//_____________________________________________________________________

Float_t CandShowerSRHandle::GetMinPhysU(Int_t iplane) const{
  Float_t minU = 999;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetPlaneView()==PlaneView::kU &&
       (ssh->GetClusterID()==ClusterType::kEMLike ||
        ssh->GetClusterID()==ClusterType::kHadLike ||
        ssh->GetClusterID()==ClusterType::kTrkLike) ) {
      Float_t testMinU = ssh->GetMinU(iplane);
      if(testMinU<minU) minU = testMinU;
    }
  }
  //if(minU==999) return 0;
  return minU;
}
//_____________________________________________________________________

Float_t CandShowerSRHandle::GetMinPhysV(Int_t iplane) const{
  Float_t minV = 999;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetPlaneView()==PlaneView::kV &&
       (ssh->GetClusterID()==ClusterType::kEMLike ||
        ssh->GetClusterID()==ClusterType::kHadLike ||
        ssh->GetClusterID()==ClusterType::kTrkLike) ) {
      Float_t testMinV = ssh->GetMinV(iplane);
      if(testMinV<minV) minV = testMinV;
    }
  }
  //if(minV==999) return 0;
  return minV;
}

//_____________________________________________________________________

Float_t CandShowerSRHandle::GetMaxPhysU(Int_t iplane) const{
  Float_t maxU = -999;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetPlaneView()==PlaneView::kU &&
       (ssh->GetClusterID()==ClusterType::kEMLike ||
        ssh->GetClusterID()==ClusterType::kHadLike ||
        ssh->GetClusterID()==ClusterType::kTrkLike) ) {
      Float_t testMaxU = ssh->GetMaxU(iplane);
      if(testMaxU>maxU) maxU = testMaxU;
    }
  }
  //if(maxU==-999) return 0;
  return maxU;
}

//_____________________________________________________________________

Float_t CandShowerSRHandle::GetMaxPhysV(Int_t iplane) const {
  Float_t maxV = -999;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetPlaneView()==PlaneView::kV &&
       (ssh->GetClusterID()==ClusterType::kEMLike ||
        ssh->GetClusterID()==ClusterType::kHadLike ||
        ssh->GetClusterID()==ClusterType::kTrkLike) ) {
      Float_t testMaxV = ssh->GetMaxV(iplane);
      if(testMaxV>maxV) maxV = testMaxV;
    }
  }
  //if(maxV==-999) return 0;
  return maxV;
}
//_____________________________________________________________________

Int_t CandShowerSRHandle::GetNPhysStrips(Int_t iplane) const{
  Int_t nstrips=0;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetClusterID()==ClusterType::kEMLike  ||
       ssh->GetClusterID()==ClusterType::kHadLike ||
       ssh->GetClusterID()==ClusterType::kTrkLike ) {
      nstrips += ssh->GetNStrips(iplane);
    }
  }
  return nstrips;
}
//_____________________________________________________________________

Bool_t CandShowerSRHandle::BelongsWithTrack(CandTrackHandle * trk, 
                                            AlgConfig & ac, 
                                            const VldContext * vldcptr, 
                                            Double_t tolTPos2,
                                            Double_t tolZPos, 
                                            Double_t tolTime){

  MSG("ShowerSR",Msg::kDebug) << "In CandShowerSRHandle::BelongsWithTrack()"
                              << endl;
  if(!trk) return false;
  Int_t pSMPlaneLast = ac.GetInt("SMPlaneLast");
  Int_t pSMPlaneFirst = ac.GetInt("SMPlaneFirst");
  VldContext vldc = *vldcptr;
  UgliGeomHandle ugh(vldc); 

  //get total charge of shower in PEs for U and V:
  Double_t shwChargeU = 0; Double_t shwChargeV = 0;
  Double_t shwStripsU = 0; Double_t shwStripsV = 0;
  Bool_t noPhysicalSubShowersU = true;
  Bool_t noPhysicalSubShowersV = true;
  for(int i=0;i<=this->GetLastSubShower();i++){    
    const CandSubShowerSRHandle *subshw = this->GetSubShower(i);
    if(subshw->GetPlaneView()==PlaneView::kU) {
      shwChargeU += subshw->GetCharge(CalStripType::kPE);
      shwStripsU += subshw->GetNStrip();
      if(subshw->GetClusterID()==ClusterType::kEMLike ||
         subshw->GetClusterID()==ClusterType::kHadLike ||
         subshw->GetClusterID()==ClusterType::kTrkLike) {
        noPhysicalSubShowersU = false;
      }
    }
    else if(subshw->GetPlaneView()==PlaneView::kV) {
      shwChargeV += subshw->GetCharge(CalStripType::kPE);
      shwStripsV += subshw->GetNStrip();
      if(subshw->GetClusterID()==ClusterType::kEMLike ||
         subshw->GetClusterID()==ClusterType::kHadLike ||
         subshw->GetClusterID()==ClusterType::kTrkLike) {
        noPhysicalSubShowersV = false;
      }
    }
  }
  Double_t shwChargeAsym = ( TMath::Abs(shwChargeU - shwChargeV) / 
                             (shwChargeU+shwChargeV) );

  //get time difference
  Double_t dt = trk->GetVtxT() - this->GetVtxT();
  if(trk->IsTPosValid(this->GetVtxPlane())){
    dt = this->GetVtxT() - trk->GetT(this->GetVtxPlane());
  }
  if(TMath::Abs(dt*1e9)>1e5) dt = 0; //can't trust it, don't use it

  //get begin and end planes
  Int_t fwdTrkPlane = min(trk->GetVtxPlane(),trk->GetEndPlane());
  Int_t bckTrkPlane = max(trk->GetVtxPlane(),trk->GetEndPlane());
  Int_t fwdShwPlane = min(this->GetVtxPlane(),this->GetEndPlane());
  Int_t bckShwPlane = max(this->GetVtxPlane(),this->GetEndPlane());

  if( bckShwPlane>=fwdTrkPlane && fwdShwPlane<=bckTrkPlane && 
      TMath::Abs(dt)<tolTime) {   // trk overlaps in Z and Time

    MSG("ShowerSR",Msg::kDebug) << " trk and shower overlap in Z and time " 
                                   << endl;

    Float_t tolTPos = 0.35; //hard coded for now

    Bool_t matchu=false;
    Bool_t matchv=false;
    Bool_t nonphysmatchu=false;
    Bool_t nonphysmatchv=false;
    // look for trk passing through shower:
    for (Int_t iplane=fwdShwPlane; iplane<=bckShwPlane;iplane++){
      PlexPlaneId plnid(GetVldContext()->GetDetector(),iplane,false);
      if(plnid.GetPlaneView()==PlaneView::kU){
        Double_t trkU = trk->GetU(iplane);
        if(!trk->IsTPosValid(iplane)){
          if(TMath::Abs(iplane - trk->GetVtxPlane()) < 
             TMath::Abs(iplane - trk->GetEndPlane())) {
            UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
            Double_t dz=scintpln.GetZ0()-trk->GetVtxZ();
            trkU = (trk->GetVtxU() + 
                    dz*trk->GetVtxDirCosU()/trk->GetVtxDirCosZ());
          }
          else{
            UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
            Double_t dz=scintpln.GetZ0()-trk->GetEndZ();
            trkU = (trk->GetEndU() + 
                    dz*trk->GetEndDirCosU()/trk->GetEndDirCosZ());
          }
        }
        if(this->GetNPhysStrips(iplane)>0 && 
           trkU>=this->GetMinPhysU(iplane)-tolTPos && 
           trkU<=this->GetMaxPhysU(iplane)+tolTPos) {
          MSG("ShowerSR", Msg::kDebug) << " u match! (details follow)" 
                                          << endl;
          matchu=true;
        }
        else if(noPhysicalSubShowersU &&
                trkU>=this->GetMinU(iplane)-tolTPos &&
                trkU<=this->GetMaxU(iplane)+tolTPos) {
          MSG("ShowerSR", Msg::kDebug) << " u match! (details follow)" 
                                          << endl;
          nonphysmatchu=true;
        }
        MSG("ShowerSR",Msg::kDebug) 
          << " plane  " << iplane << " trk u " << trkU 
          << " min/max U " << this->GetMinU(iplane) << "/" 
          << this->GetMaxU(iplane)
          << " phys min/max U " << this->GetMinPhysU(iplane) << "/" 
          << this->GetMaxPhysU(iplane) << endl;
      }
      else if(plnid.GetPlaneView()==PlaneView::kV){
        Double_t trkV = trk->GetV(iplane);
        if(!trk->IsTPosValid(iplane)){
          if(TMath::Abs(iplane - trk->GetVtxPlane()) < 
             TMath::Abs(iplane - trk->GetEndPlane())) {
            UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
            Double_t dz=scintpln.GetZ0()-trk->GetVtxZ();
            trkV = (trk->GetVtxV() + 
                    dz*trk->GetDirCosV()/trk->GetDirCosZ());
          }
          else{
            UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
            Double_t dz=scintpln.GetZ0()-trk->GetEndZ();
            trkV= (trk->GetEndV() + 
                   dz*trk->GetEndDirCosV()/trk->GetEndDirCosZ());
          }
        }
        if(this->GetNPhysStrips(iplane)>0 &&
           trkV>=this->GetMinPhysV(iplane)-tolTPos && 
           trkV<=this->GetMaxPhysV(iplane)+tolTPos) {
          MSG("ShowerSR", Msg::kDebug) << " v match! (details follow)" 
                                          << endl;
          matchv=true;
        }
        else if(noPhysicalSubShowersV &&
                trkV>=this->GetMinV(iplane)-tolTPos &&
                trkV<=this->GetMaxV(iplane)+tolTPos) {
          MSG("ShowerSR", Msg::kDebug) << " v match! (details follow)" 
                                          << endl;
          nonphysmatchv=true;
        }
        MSG("ShowerSR",Msg::kDebug) 
          << " plane  " << iplane << " trk v " << trkV 
          << " min/max V " << this->GetMinV(iplane) << "/" 
          << this->GetMaxV(iplane)
          << " phys min/max V " << this->GetMinPhysV(iplane) << "/" 
          << this->GetMaxPhysV(iplane) << endl;
      }
        
      if(matchu && matchv) {
        MSG("ShowerSR", Msg::kDebug) << " Phys 3D  match! " << endl;
        return true;
      }
      else if( matchu || matchv ){ //look for small asymmetric shws along trk
        MSG("ShowerSR", Msg::kDebug) 
          << "Match Stat: " << matchu << " " << matchv << " "
          << "Shower Charge: " << shwChargeU << " " << shwChargeV << " "
          << "Shower Charge Asymmetry: " << shwChargeAsym << endl;
        if(matchu && shwChargeV<10. && 
           shwChargeV/shwStripsV < 3. && shwChargeAsym>0.7 ) return true;
        else if(matchv && shwChargeU<10. && 
                shwChargeU/shwStripsU < 3. && shwChargeAsym>0.7 ) return true;
        else if( (matchu && nonphysmatchv) || 
                 (matchv && nonphysmatchu) ) return true;
      }
      else if(nonphysmatchu && nonphysmatchv) {
        MSG("ShowerSR", Msg::kDebug) << " Non-Phys 3D  match! " << endl;
        return true;
      }
      else if(nonphysmatchu || nonphysmatchv) {
        MSG("ShowerSR", Msg::kDebug) 
          << "Match Stat: " << nonphysmatchu << " " << nonphysmatchv << " "
          << "Shower Charge: " << shwChargeU << " " << shwChargeV << " "
          << "Shower Charge Asymmetry: " << shwChargeAsym << endl;
        if(shwChargeU/shwStripsU < 3. || shwChargeV/shwStripsV < 3. ||
           shwChargeAsym > 0.6) return true;
      }
    }

    //loop over all strips in shower
    //calculate ph weighted distance from track
    Double_t phDistance = 0;
    Double_t phDistanceU = 0;
    Double_t phDistanceV = 0;
    TIter stripItr(GetDaughterIterator());
    while (CandStripHandle *strip = 
           dynamic_cast<CandStripHandle*>(stripItr())) {
      if(strip->GetPlaneView()==PlaneView::kU) {
        if(!trk->IsTPosValid(strip->GetPlane())) continue;
        Double_t trkU = trk->GetU(strip->GetPlane());
        phDistance += strip->GetCharge(CalDigitType::kPE) * 
          TMath::Power(strip->GetTPos()-trkU,2);
        phDistanceU += strip->GetCharge(CalDigitType::kPE) * 
          TMath::Power(strip->GetTPos()-trkU,2);
      }
      else if(strip->GetPlaneView()==PlaneView::kV){
        if(!trk->IsTPosValid(strip->GetPlane())) continue;
        Double_t trkV = trk->GetV(strip->GetPlane());
        phDistance += strip->GetCharge(CalDigitType::kPE) * 
          TMath::Power(strip->GetTPos()-trkV,2);
        phDistanceV += strip->GetCharge(CalDigitType::kPE) * 
          TMath::Power(strip->GetTPos()-trkV,2);
      }
    }
    if(shwChargeU+shwChargeV>0) phDistance /= (shwChargeU+shwChargeV);
    phDistance = TMath::Sqrt(phDistance);
    if(shwChargeU>0) phDistanceU /= shwChargeU;
    phDistanceU = TMath::Sqrt(phDistanceU);
    if(shwChargeV>0) phDistanceV /= shwChargeV;
    phDistanceV = TMath::Sqrt(phDistanceV);
    MSG("ShowerSR",Msg::kDebug) << "phDistance = " << phDistance 
                                   << " phDistanceU = " << phDistanceU
                                   << " phDistanceV = " << phDistanceV
                                   << endl;
    if(phDistance  < tolTPos  || 
       (phDistanceU < tolTPos && shwChargeAsym>0.6) || 
       (phDistanceV < tolTPos && shwChargeAsym>0.6) ) return true;
  }

  //strip based checks failed!
  //Try comparing shower vertex with track trajectory at shw vertex plane
  if( this->GetVtxZ()>min(trk->GetVtxZ(),trk->GetEndZ())-tolZPos  &&
      this->GetVtxZ()<max(trk->GetVtxZ(),trk->GetEndZ())+tolZPos) {
    Double_t du = trk->GetVtxU() - this->GetVtxU();
    Double_t dv = trk->GetVtxV() - this->GetVtxV();     
    dt = trk->GetVtxT() - this->GetVtxT();
    if(trk->IsTPosValid(this->GetVtxPlane())){
      du = this->GetVtxU()-trk->GetU(this->GetVtxPlane());
      dv = this->GetVtxV()-trk->GetV(this->GetVtxPlane());
      dt = this->GetVtxT()-trk->GetT(this->GetVtxPlane());
    }
    MSG("ShowerSR",Msg::kDebug) 
      << "Vertex check at plane " << this->GetVtxPlane() 
      << ": du^{2}+dv^{2} = " << du*du+dv*dv << "(/" << tolTPos2 
      << "); dt = " << dt*1.e9 << "(/" << tolTime*1e9 << ")" << endl;
    // check whether vertex within minimum distance from trk 
    if(dt<tolTime && du*du+dv*dv<tolTPos2) return true;   
  }
  
  // finally, check distance between vertices 
  Double_t trkZ = trk->GetVtxZ();
  Double_t trkU = trk->GetVtxU();
  Double_t trkV = trk->GetVtxV();
  Double_t trkT = trk->GetVtxT();
  Int_t trkPlane = trk->GetVtxPlane();
  Double_t trkDirU = trk->GetVtxDirCosU();
  Double_t trkDirV = trk->GetVtxDirCosV();
  Double_t trkDirZ = trk->GetVtxDirCosZ();

  if(TMath::Abs(this->GetVtxZ() - trk->GetEndZ()) < 
     TMath::Abs(this->GetVtxZ() - trk->GetVtxZ())){
    trkZ = trk->GetEndZ();
    trkU = trk->GetEndU();
    trkV = trk->GetEndV();
    trkT = trk->GetEndT();
    trkDirU = trk->GetEndDirCosU();
    trkDirV = trk->GetEndDirCosV();
    trkDirZ = trk->GetEndDirCosZ();
    trkPlane = trk->GetEndPlane();
  }

  Double_t dz = trkZ - this->GetVtxZ();
  Double_t du = trkU - this->GetVtxU();
  if(trkDirZ>0) du = (trkU - dz*trkDirU/trkDirZ) - this->GetVtxU();
  Double_t dv = trkV - this->GetVtxV();
  if(trkDirZ>0) dv = (trkV - dz*trkDirV/trkDirZ) - this->GetVtxV();
  dt = trkT - this->GetVtxT(); 

  //compensate dz for SM gap if necessary
  if (vldcptr->GetDetector()==Detector::kFar) {
    if((trkPlane<=pSMPlaneLast && this->GetVtxPlane()>=pSMPlaneFirst) ||
       (trkPlane>=pSMPlaneLast && this->GetVtxPlane()<=pSMPlaneFirst)) 
      dz = (trkPlane - this->GetVtxPlane())*0.0594;
  }

  MSG("ShowerSR",Msg::kDebug) 
    << "Delta(vertex) shower/track (u,v,z,t): " 
    << du << " " << dv << " " << dz <<" "  << dt*1.e9 << endl;
  if (du*du+dv*dv<tolTPos2 && TMath::Abs(dz)<tolZPos && 
      TMath::Abs(dt)<tolTime) return true;
  
  // now try shower End
  dz = trkZ - this->GetEndZ();
  du = (trkU - dz*trkDirU/trkDirZ) - this->GetEndU();
  dv = (trkV - dz*trkDirV/trkDirZ) - this->GetEndV();
  dt = trkT - this->GetVtxT(); 

  //compensate dz for SM gap if necessary
  if (vldcptr->GetDetector()==Detector::kFar) {
    if((trkPlane<=pSMPlaneLast && this->GetEndPlane()>=pSMPlaneFirst) ||
       (trkPlane>=pSMPlaneLast && this->GetEndPlane()<=pSMPlaneFirst)) 
      dz = (trkPlane - this->GetEndPlane())*0.0594;
  }

  MSG("ShowerSR",Msg::kDebug) 
    << "Delta(vertex) shower/track (u,v,z,t): " 
    << du << " " << dv << " " << dz <<" "  << dt*1.e9 << endl;
  if (du*du+dv*dv<tolTPos2 && TMath::Abs(dz)<tolZPos && 
      TMath::Abs(dt)<tolTime) return true;

  return false;
}
//_____________________________________________________________________

Bool_t CandShowerSRHandle::BelongsWithShower(CandShowerHandle * shw, 
                                             AlgConfig & /*ac*/, 
                                             const VldContext * /*vldcptr*/, 
                                             Double_t tolTPos2,
                                             Double_t tolZPos, 
                                             Double_t tolTime){

  MSG("ShowerSR",Msg::kDebug) << "In BelongsWithShower()" << endl;

  tolZPos  = 10.0; //(m)
  tolTime  = 7.5; //(ns)
  tolTPos2 = 10.0*10.0; //(m^2)

  if(shw->InheritsFrom("CandShowerSRHandle")) {
    CandShowerSRHandle *shwSR = dynamic_cast<CandShowerSRHandle*> (shw);
    //first try to merge showers very nearby in space and time:
    Double_t xyzDistance = 
      TMath::Sqrt( ( TMath::Power(shwSR->GetVtxU()-this->GetVtxU(),2) +
                     TMath::Power(shwSR->GetVtxV()-this->GetVtxV(),2) +
                     TMath::Power(shwSR->GetVtxZ()-this->GetVtxZ(),2) ) );
    Double_t temp = 0;
    Double_t xyztDistance =
      TMath::Sqrt( ( TMath::Power(shwSR->GetVtxU()-this->GetVtxU(),2) +
                     TMath::Power(shwSR->GetVtxV()-this->GetVtxV(),2) +
                     TMath::Power(shwSR->GetVtxZ()-this->GetVtxZ(),2) +
                     TMath::Power(Mphysical::c_light*
                                  (shwSR->GetPHWTime(temp) - 
                                   this->GetPHWTime(temp)),2) ) );
    
    MSG("ShowerSR",Msg::kDebug) << "xyzDistance = " << xyzDistance << endl;
    MSG("ShowerSR",Msg::kDebug) << "xyztDistance = " << xyztDistance << endl;
    
    if(xyztDistance<1.5) return true;
  }

  Double_t minZDistance = shw->GetVtxZ() - this->GetVtxZ();
  if(minZDistance<0) minZDistance = this->GetVtxZ() - shw->GetVtxZ();
  Double_t minTDistance2 = (TMath::Power(shw->GetVtxU() - this->GetVtxU(),2) +
                            TMath::Power(shw->GetVtxV() - this->GetVtxV(),2) );

  if( minZDistance > tolZPos ) return false;
  if( minTDistance2 > tolTPos2 ) return false;

  MSG("ShowerSR",Msg::kDebug) << "Z Separation = " << minZDistance
                              << " (< " << tolZPos << ")" << endl;
  MSG("ShowerSR",Msg::kDebug) << "TPos^2 Separation = " << minTDistance2
                                 << " (< " << tolTPos2 << ")" << endl;

  Double_t realDistance = TMath::Sqrt(TMath::Power(minZDistance,2) + 
                                      minTDistance2 );
  MSG("ShowerSR",Msg::kDebug) << "True Separation = " 
                              << realDistance << endl;

  //assume upper limit to be 1-sigma less energy reco'd than was really there
  Double_t uppShwEn = TMath::Power(0.55/2. + 
           TMath::Sqrt(0.55*0.55/4. + shw->GetEnergy(CandShowerHandle::kCC)),2);

  Double_t mOverE = (Mphysical::neutron_mass_c2/Munits::hep2baseMass)/ 
    shw->GetEnergy(CandShowerHandle::kCC);
  Double_t upp_mOverE = (Mphysical::neutron_mass_c2/Munits::hep2baseMass)/uppShwEn;

  if(mOverE > 0.98) mOverE = upp_mOverE;
  if(mOverE > 0.98) mOverE = 0.98;  // <=> ~200 MeV neutron  
  Double_t neutronSpeed = ( TMath::Sqrt(1. - TMath::Power(mOverE,2)) * 
                            Mphysical::c_light );

  Double_t predTime = 1.0e9*realDistance/neutronSpeed;
  Double_t predTime_c = 1.0e9*realDistance/Mphysical::c_light;
  if(shw->GetVtxZ() - this->GetVtxZ()<0) {
    predTime = -predTime; //if shw is behind this
    predTime_c = -predTime_c;
  }
  
  Double_t realPHWTime = 0;
  Double_t realAveTime = 0;
  Double_t phErr1=0,phErr2=0,err1=0,err2=0;
  if(shw->InheritsFrom("CandShowerSRHandle")) {
    CandShowerSRHandle *shwSR = dynamic_cast<CandShowerSRHandle*> (shw);
    realPHWTime = 1.0e9*(shwSR->GetPHWTime(phErr2) - this->GetPHWTime(phErr1));
    realAveTime = 1.0e9*(shwSR->GetAveTime(err2) - this->GetAveTime(err1));
    phErr1*=1e9; phErr2*=1e9; err1*=1e9; err2*=1e9;
  }
  else {
    realPHWTime = 1.0e9*(shw->GetVtxT() - this->GetEndT());
    realAveTime = realPHWTime;
  }

  Double_t realPHWTimeErr = TMath::Sqrt(phErr1*phErr1 + phErr2*phErr2);
  Double_t realAveTimeErr = TMath::Sqrt(err1*err1 + err2*err2);

  MSG("ShowerSR",Msg::kDebug) 
    << "\nPredicted Time Separation = " << predTime 
    << " (assuming c = " << predTime_c << ")"
    << "\nActual Time Separation = (PH-weighted average,simple average) " 
    << realPHWTime << " +/- " << realPHWTimeErr
    << " " << realAveTime << " +/- " << realAveTimeErr << endl;
  
  if(TMath::Abs(predTime-realPHWTime)<realPHWTimeErr ||
     TMath::Abs(predTime-realAveTime)<realAveTimeErr ||
     ( realPHWTime+realPHWTimeErr>=predTime_c && 
       realPHWTime-realPHWTimeErr<=predTime) || 
     ( realAveTime+realAveTimeErr>=predTime_c && 
       realAveTime-realAveTimeErr<=predTime) ) return true;
  return false;
  
}

//_____________________________________________________________________

Bool_t CandShowerSRHandle::BuriedTrack(CandTrackHandle *track,
                                       Double_t pMinTrackIsolation,
                                       Int_t *stats)
{
  if(stats) {
    for(int i=0;i<8;i++) stats[i] = -1;
  }

  Int_t firstTrkPlane = TMath::Min(track->GetBegPlane(),
                                   track->GetEndPlane());
  Int_t lastTrkPlane = TMath::Max(track->GetBegPlane(),
                                  track->GetEndPlane());
  Int_t firstShwPlane = TMath::Min(this->GetBegPlane(),
                                   this->GetEndPlane());
  Int_t lastShwPlane = TMath::Max(this->GetBegPlane(),
                                  this->GetEndPlane());

  Int_t firstSharedPlane = 500;
  Int_t lastSharedPlane = 0;

  CandFitTrackHandle *fittrack = NULL;
  if(track->InheritsFrom("CandFitTrackHandle")) {
    fittrack = dynamic_cast<CandFitTrackHandle*> (track);    
    MSG("ShowerSR",Msg::kDebug) << "Track Pass? " << fittrack->GetPass() << endl;
  }

  //look at number of track planes protruding from shower plane region
  Int_t nSharedPlanes = lastShwPlane - firstTrkPlane + 1;
  if(nSharedPlanes<0) nSharedPlanes = 0;
  else {
    if(firstTrkPlane<firstShwPlane) {
      nSharedPlanes += firstTrkPlane-firstShwPlane;
      firstSharedPlane = firstShwPlane;
    }
    else firstSharedPlane = firstTrkPlane;
    if(lastTrkPlane<lastShwPlane) {
      nSharedPlanes += lastTrkPlane-lastShwPlane;
      lastSharedPlane = lastTrkPlane;
    }
    else lastSharedPlane = lastShwPlane;
  }
  if(nSharedPlanes<=0) return false;
  if(stats) stats[0] = nSharedPlanes;

  MSG("ShowerSR",Msg::kDebug) << "lastTrkPlane = " << lastTrkPlane 
                              << " firstTrkPlane = " << firstTrkPlane     
                              << endl;
  MSG("ShowerSR",Msg::kDebug) << "lastShwPlane = " << lastShwPlane
                              << " firstShwPlane = " << firstShwPlane 
                              << endl;
  MSG("ShowerSR",Msg::kDebug) << "nSharedPlanes = " << nSharedPlanes << endl;

  Int_t nMissingIsoPlanes = 0;
  Int_t nMissingSharedPlanes = 0;
  Float_t planeCharge[500] = {};
  if(fittrack) {
    if(fittrack->GetPass()) {
      for(int iplane=firstTrkPlane;iplane<=lastTrkPlane;iplane++){
        MSG("ShowerSR",Msg::kVerbose) << iplane << " " 
                                      << track->GetPlaneCharge(iplane) << endl;
        planeCharge[iplane] = track->GetPlaneCharge(iplane);
        if(track->GetPlaneCharge(iplane)==0) {
          if(iplane>=firstSharedPlane && iplane<=lastSharedPlane) 
            nMissingSharedPlanes+=1;
          else nMissingIsoPlanes+=1;
        }
      }
    }
    else {
      const CandSliceHandle *slice = fittrack->GetCandSlice();
      nMissingSharedPlanes = lastSharedPlane - firstSharedPlane + 1;
      nMissingIsoPlanes = lastTrkPlane - firstTrkPlane + 1 - nMissingSharedPlanes;
      for(int iplane=firstTrkPlane;iplane<=lastTrkPlane;iplane++){
        TIter stripitr(slice->GetDaughterIterator());
        while(CandStripHandle *strip =
              dynamic_cast<CandStripHandle*>(stripitr())){
          if(strip->GetPlane()==iplane) {
            PlexPlaneId plnid(track->GetVldContext()->GetDetector(),
                              iplane,false);
            if(plnid.GetPlaneView() == PlaneView::kU){
              Float_t trkU = track->GetU(iplane);
              if(trkU<strip->GetTPos()+0.0205 && trkU>strip->GetTPos()-0.0205) {
                planeCharge[iplane] = strip->GetCharge();
                if(iplane>=firstSharedPlane && iplane<=lastSharedPlane)
                  nMissingSharedPlanes-=1;
                else nMissingIsoPlanes-=1;
              }
            }
            else if(plnid.GetPlaneView() == PlaneView::kV){
              Float_t trkV = track->GetV(iplane);
              if(trkV<strip->GetTPos()+0.0205 && trkV>strip->GetTPos()-0.0205) {
                planeCharge[iplane] = strip->GetCharge();
                if(iplane>=firstSharedPlane && iplane<=lastSharedPlane)
                  nMissingSharedPlanes-=1;
                else nMissingIsoPlanes-=1;
              }
            }
          }
        }
      }
    }
  }

  Int_t nIsolatedPlanes = ( lastTrkPlane - firstTrkPlane + 1 -
                            nMissingIsoPlanes - nSharedPlanes);

  MSG("ShowerSR",Msg::kDebug) << "nMissingIsoPlanes = " 
                              << nMissingIsoPlanes
                              << " nMissingSharedPlanes = " 
                              << nMissingSharedPlanes 
                              << " nIsolatedPlanes = " 
                              << nIsolatedPlanes << endl;

  if(stats) {
    stats[1] = nMissingIsoPlanes;
    stats[2] = nMissingSharedPlanes;
  }
  if(nIsolatedPlanes>=pMinTrackIsolation) return false;

  //now look transversely at overlapping planes:
  Int_t nBuriedPlanes = 0;
  Int_t nPhysBuriedPlanes = 0;
  for(int iplane=firstSharedPlane;iplane<=lastSharedPlane;iplane++){
    PlexPlaneId plnid(track->GetVldContext()->GetDetector(),iplane,false);
    if(planeCharge[iplane]==0) continue;
    if(plnid.GetPlaneView() == PlaneView::kU){
      Float_t minShw = this->GetMinU(iplane)-0.02;
      Float_t maxShw = this->GetMaxU(iplane)+0.02;
      Float_t minPhysShw = this->GetMinPhysU(iplane)-0.02;
      Float_t maxPhysShw = this->GetMaxPhysU(iplane)+0.02;
      Float_t trkU = track->GetU(iplane);
      if(trkU>minShw&&trkU<maxShw) nBuriedPlanes+=1;
      if(this->GetNPhysStrips(iplane)&&
         trkU>minPhysShw&&trkU<maxPhysShw) nPhysBuriedPlanes+=1;
    }
    else if(plnid.GetPlaneView() == PlaneView::kV){
      Float_t minShw = this->GetMinV(iplane)-0.02;
      Float_t maxShw = this->GetMaxV(iplane)+0.02;
      Float_t minPhysShw = this->GetMinPhysV(iplane)-0.02;
      Float_t maxPhysShw = this->GetMaxPhysV(iplane)+0.02;
      Float_t trkV = track->GetV(iplane);
      if(trkV>minShw&&trkV<maxShw) nBuriedPlanes+=1;
      if(this->GetNPhysStrips(iplane)&&
         trkV>minPhysShw&&trkV<maxPhysShw) nPhysBuriedPlanes+=1;
    }
  }

  nIsolatedPlanes = ( lastTrkPlane - firstTrkPlane + 1 - 
                      nMissingIsoPlanes - nBuriedPlanes - 
                      nMissingSharedPlanes);
  MSG("ShowerSR",Msg::kDebug) << "nBuriedPlanes = " << nBuriedPlanes 
                              << " nIsolatedPlanes = " << nIsolatedPlanes << endl;
  if(stats) {
    stats[3] = nBuriedPlanes;
    stats[4] = nPhysBuriedPlanes;
  }

  if(nIsolatedPlanes>=pMinTrackIsolation) return false;

  //count up number of times track passes 
  //through x-talk/low ph halo hits and track-like hits
  Int_t nXTalk = 0;
  Int_t nTrkLike = 0;
  Int_t nTrkLikePH = 0;  
  for(int iplane=firstSharedPlane;iplane<=lastSharedPlane;iplane++){
    PlexPlaneId plnid(track->GetVldContext()->GetDetector(),iplane,false);  
    for(int i=0;i<=this->GetLastSubShower();i++){
      const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
      if(ssh->GetPlaneView()!=plnid.GetPlaneView()) continue;
      if(iplane<ssh->GetBegPlane()||iplane>ssh->GetEndPlane()) continue;
      TIter stripItr(ssh->GetDaughterIterator());
      while(CandStripHandle *strip = 
            dynamic_cast<CandStripHandle*>(stripItr())) {
        bool match = false;
        if(strip->GetPlane()!=iplane) continue;
        if(plnid.GetPlaneView()==PlaneView::kU && 
           TMath::Abs(strip->GetTPos()-track->GetU(iplane))<0.0205){
          match = true;
        }
        else if(plnid.GetPlaneView()==PlaneView::kV && 
                TMath::Abs(strip->GetTPos()-track->GetV(iplane))<0.0205){
          match = true;
        }
        if(match){
          if(ssh->GetClusterID()==ClusterType::kXTalk) nXTalk+=1;
          else if(ssh->GetClusterID()==ClusterType::kHalo && 
                  strip->GetCharge(CalDigitType::kPE)<2.0) nXTalk+=1;
          else if(ssh->GetClusterID()==ClusterType::kTrkLike) nTrkLike+=1;        
          else if(ssh->GetClusterID()==ClusterType::kEMLike ||
                  ssh->GetClusterID()==ClusterType::kHadLike) {
            TIter stripItr2(this->GetDaughterIterator());
            Bool_t isolatedTrackHit = true;
            while(CandStripHandle *strip2 = 
                  dynamic_cast<CandStripHandle*>(stripItr2())) {
              if(strip2->GetPlane()!=iplane) continue;        
              if(plnid.GetPlaneView()==PlaneView::kU) {
                if( TMath::Abs(strip2->GetTPos()-track->GetU(iplane)-0.041)<0.0205 ||
                    TMath::Abs(strip2->GetTPos()-track->GetU(iplane)+0.041)<0.0205 )
                  isolatedTrackHit = false;
              }
              else if(plnid.GetPlaneView()==PlaneView::kV){
                if( TMath::Abs(strip2->GetTPos()-track->GetV(iplane)-0.041)<0.0205 ||
                    TMath::Abs(strip2->GetTPos()-track->GetV(iplane)+0.041)<0.0205 )
                  isolatedTrackHit = false;
              }
            }
            if(isolatedTrackHit) nTrkLikePH++;      
          }
          break;
        }
      }
    }
  }

  //add back in trk-like planes:
  nIsolatedPlanes += nTrkLike;
  MSG("ShowerSR",Msg::kDebug) << "nTrkLike = " << nTrkLike
                              << " nTrkLikePH = " << nTrkLikePH
                              << " nIsolatedPlanes = " << nIsolatedPlanes << endl;
  if(stats) {
    stats[5] = nTrkLike;
    stats[6] = nTrkLikePH;
    stats[7] = nXTalk;
  }
  if(nIsolatedPlanes+nTrkLikePH>=pMinTrackIsolation) return false;

  //now consider only physics subshowers to decide if track
  //is deep within shower core, or if it's in the halo of the shower
  //don't allow xtalk level hits in the halo to be included in the count
  nIsolatedPlanes += (nBuriedPlanes - nPhysBuriedPlanes - nXTalk);
  MSG("ShowerSR",Msg::kDebug) << "nPhysBuriedPlanes = " << nPhysBuriedPlanes 
                              << " nIsolatedPlanes = " << nIsolatedPlanes << endl;
  if(nIsolatedPlanes>=pMinTrackIsolation) return false;
  
  return true;

}

//_____________________________________________________________________

Bool_t CandShowerSRHandle::IsUnphysical(Float_t xtalkFrac,
                                        Float_t xtalkCut)
{
  Int_t nPhysSS = 0;
  for(int i=0;i<=this->GetLastSubShower();i++){
    const CandSubShowerSRHandle *ssh = this->GetSubShower(i);
    if(ssh->GetClusterID()==ClusterType::kEMLike ||
       ssh->GetClusterID()==ClusterType::kHadLike ||
       ssh->GetClusterID()==ClusterType::kTrkLike) nPhysSS+=1;
  }
  if(nPhysSS>1) return false;
  Int_t nXTlk = 0;
  Int_t nSig = 0;
  Float_t aveCharge = 0;
  Float_t maxTPosU = -4;
  Float_t minTPosU = +4;
  Float_t maxTPosV = -4;
  Float_t minTPosV = +4;
  TIter stripItr(this->GetDaughterIterator());
  while (CandStripHandle *strip = 
         dynamic_cast<CandStripHandle*>(stripItr())) {  
    aveCharge += strip->GetCharge(CalDigitType::kPE);
    if(strip->GetCharge(CalDigitType::kPE)<xtalkCut) nXTlk+=1;
    if(strip->GetCharge(CalDigitType::kPE)>5.0) nSig+=1;
    if(strip->GetPlaneView()==PlaneView::kU) {
      if(strip->GetTPos()>maxTPosU) maxTPosU = strip->GetTPos();
      if(strip->GetTPos()<minTPosU) minTPosU = strip->GetTPos();
    }
    else if(strip->GetPlaneView()==PlaneView::kV){
      if(strip->GetTPos()>maxTPosV) maxTPosV = strip->GetTPos();
      if(strip->GetTPos()<minTPosV) minTPosV = strip->GetTPos();
    }
  }
  if(nSig>1) return false;
  if(Float_t(nXTlk)/Float_t(this->GetNStrip())>xtalkFrac) return true;
  else if(aveCharge/Float_t(this->GetNStrip())<xtalkCut) return true;
  return false;
}

//_____________________________________________________________________

Double_t CandShowerSRHandle::GetPHWTime(Double_t &error)
{
  //find the pulse height weighted mean and rms of strip time
  //then refine estimate by removing outliers > 1 rms from mean
  Double_t PHWTime = 0.;
  Double_t PHTot = 0.;
  Double_t AveTime = 0.;
  Double_t TimeSq = 0.;
  Double_t Tot = 0.;
  TIter stripItr(GetDaughterIterator());
  while (CandStripHandle *strip = 
         dynamic_cast<CandStripHandle*>(stripItr())) {  
    Tot += 1.;
    AveTime += strip->GetTime();
    TimeSq += strip->GetTime()*strip->GetTime();
  }
  Double_t mean = AveTime/Tot;
  Double_t rms = 0;
  if(TimeSq/Tot - mean*mean>0) rms = TMath::Sqrt(TimeSq/Tot - mean*mean);
  TimeSq = 0;
  stripItr.Reset();
  while (CandStripHandle *strip = 
         dynamic_cast<CandStripHandle*>(stripItr())) {
    if(TMath::Abs(strip->GetTime() - mean) < rms) {
      PHTot += strip->GetCharge(CalDigitType::kPE);
      PHWTime += strip->GetCharge(CalDigitType::kPE) * 
        strip->GetTime();
      TimeSq += strip->GetCharge(CalDigitType::kPE) * 
        strip->GetTime()*strip->GetTime();
    }
  }
  error = -1;
  if(PHTot<=0) return 0;
  mean = PHWTime/PHTot;
  if(TimeSq/PHTot - mean*mean>0) 
    error = TMath::Sqrt(TimeSq/PHTot - mean*mean);
  return mean;
}
//_____________________________________________________________________

Double_t CandShowerSRHandle::GetAveTime(Double_t &error)
{
  //find the mean and rms of strip time
  //then refine estimate by removing outliers > 1 rms from mean
  Double_t AveTime = 0.;
  Double_t TimeSq = 0.;
  Double_t Tot = 0.;
  TIter stripItr(GetDaughterIterator());
  while (CandStripHandle *strip = 
         dynamic_cast<CandStripHandle*>(stripItr())) {
    Tot += 1.;
    AveTime += strip->GetTime();
    TimeSq += strip->GetTime()*strip->GetTime();
  }
  if(Tot<=0) return 0;
  stripItr.Reset();
  Double_t mean = AveTime/Tot;
  Double_t rms = 0;
  if(TimeSq/Tot - mean*mean>0) rms = TMath::Sqrt(TimeSq/Tot - mean*mean);
  TimeSq = 0;
  AveTime = 0;
  Tot = 0;
  while (CandStripHandle *strip = 
         dynamic_cast<CandStripHandle*>(stripItr())) {  
    if(TMath::Abs(strip->GetTime()-mean) < rms){
      Tot += 1.;
      AveTime += strip->GetTime();
      TimeSq += strip->GetTime()*strip->GetTime();
    }
  }
  error = -1;
  if(Tot<=0) return 0;
  mean = AveTime/Tot;
  if(TimeSq/Tot - mean*mean>0) 
    error = TMath::Sqrt(TimeSq/Tot - mean*mean);
  return mean;
}

XXXITRIMP(CandShowerSRHandle)

---