AlgEventSSList.cxx

Go to the documentation of this file.

// $Id: AlgEventSSList.cxx,v 1.16 2009/02/28 22:00:32 gmieg Exp $
//
// AlgEventSSList.cxx
//
// AlgEventSSList is a concrete EventSSList Algorithm class.
//
// Author:  R. Lee 2001.03.07

#include <cassert>
#include <cmath>
#include <map>

#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
#include "CandEventSR/AlgEventSSList.h"
#include "RecoBase/CandEvent.h"
#include "RecoBase/CandEventHandle.h"
#include "RecoBase/CandEventList.h"
#include "RecoBase/CandEventListHandle.h"
#include "Conventions/PlaneView.h"
#include "CandShowerSR/CandShowerSR.h"
#include "CandShowerSR/CandShowerSRHandle.h"
#include "CandShowerSR/CandShowerSRListHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandShower.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Plex/PlexPlaneId.h"
#include "RecoBase/CandEventHandle.h"
#include "RecoBase/CandRecoHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "CandSubShowerSR/CandSubShowerSR.h"
#include "CandSubShowerSR/CandSubShowerSRHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "CandSubShowerSR/CandSubShowerSRListHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "Validity/VldContext.h"
#include "CandFitTrackSR/CandFitTrackSRHandle.h"
#include "Calibrator/Calibrator.h"
#include "CandChop/ChopHelper.h"
#include "CandChop/ChopHelp.h"
ClassImp(AlgEventSSList)

#define PITCHINMETRES 0.0594

//______________________________________________________________________
CVSID("$Id: AlgEventSSList.cxx,v 1.16 2009/02/28 22:00:32 gmieg Exp $");

//______________________________________________________________________
AlgEventSSList::AlgEventSSList()
{
}

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

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

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

  CandContext cxx(this,cx.GetMom());
  AlgFactory &af = AlgFactory::GetInstance();
  const char *pEventAlgConfig = 0;
  ac.Get("EventAlgConfig",pEventAlgConfig);
  AlgHandle ah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
  // grab all needed algorithm parameters  

  //generic parameters:
  Double_t pTrkTrkDtpos2 = ac.GetDouble("TrkTrkDtpos2");
  Double_t pTrkTrkDz = ac.GetDouble("TrkTrkDz");
  Double_t pTrkTrkDt = ac.GetDouble("TrkTrkDt");
  Double_t pShwTrkDtpos2 = ac.GetDouble("ShwTrkDtpos2");
  Double_t pShwTrkDz = ac.GetDouble("ShwTrkDz");
  Double_t pShwTrkDt = ac.GetDouble("ShwTrkDt");
  Double_t pShwShwDtpos2 = ac.GetDouble("ShwShwDtpos2");
  Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
  Double_t pShwShwDt = ac.GetDouble("ShwShwDt");
  Double_t pMaxNewShwLen = ac.GetDouble("MaxNewShwLen");
  Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
  Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
  Double_t pHitAssocMaxDist2 = ac.GetDouble("HitAssocMaxDist2");
  Double_t pMinTrackIsolation = ac.GetDouble("MinTrackIsolation");

  //SS specific parameters
  Double_t pHardBuriedFrac = ac.GetDouble("HardBuriedFrac");
  Double_t pSoftBuriedFrac = ac.GetDouble("SoftBuriedFrac");
  Double_t pDCosUVCut = ac.GetDouble("DCosUVCut");
  Double_t pTrkGapFracCut = ac.GetDouble("TrkGapFracCut");
  Double_t pTrkAsymCut = ac.GetDouble("TrkAsymCut");
  Double_t pTrkXTalkFracCut = ac.GetDouble("TrkXTalkFracCut");
  Double_t pTrkXTalkDef = ac.GetDouble("TrkXTalkDef");
  Double_t pShwXTalkFracCut = ac.GetDouble("ShwXTalkFracCut");
  Double_t pShwXTalkDef = ac.GetDouble("ShwXTalkDef");
  Int_t useChopper = ac.GetInt("UseChopper");

  MSG("EventSS",Msg::kDebug) << pHardBuriedFrac << " " << pSoftBuriedFrac << " "
                            << pDCosUVCut << " " << pTrkGapFracCut << " "
                            << pTrkAsymCut << " " << pTrkXTalkFracCut << " "
                            << pTrkXTalkDef << " " << pShwXTalkFracCut << " "
                            << pShwXTalkDef << endl;

  // grab all needed input lists
  const CandRecord *candrec = cx.GetCandRecord();
  assert(candrec);
  const VldContext *vldcptr = candrec->GetVldContext();
  assert(vldcptr);
  VldContext vldc = *vldcptr;
  const CandSliceListHandle *slicelist = 0;
  CandShowerListHandle *showerlist = 0;
  const CandTrackListHandle *tracklist = 0;
  CandSubShowerSRListHandle *subshowerlist = 0;
  const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
  for (Int_t i=0; i<=cxin->GetLast(); i++) {
    TObject *tobj = cxin->At(i);
    if (tobj->InheritsFrom("CandSliceListHandle")) {
      slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
    }
    if (tobj->InheritsFrom("CandShowerListHandle")) {
      showerlist = dynamic_cast<CandShowerListHandle*>(tobj);
    }
    if (tobj->InheritsFrom("CandTrackListHandle")) {
      tracklist = dynamic_cast<CandTrackListHandle*>(tobj);
    }
    if (tobj->InheritsFrom("CandSubShowerSRListHandle")) {
      subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj);
    }
  }
  if (!slicelist) {
    MSG("EventSS",Msg::kError) << "CandSliceListHandle missing\n";
    return;
  }

  CandShowerHandleItr * showerItr=0;
  // create iterators for track, shower, and cluster lists, keying on slice
  if (showerlist && showerlist->GetNDaughters()>0) {
    showerItr = new CandShowerHandleItr(showerlist->GetDaughterIterator());
    CandShowerHandleKeyFunc *showerKf = showerItr->CreateKeyFunc();
    showerKf->SetFun(CandShowerHandle::KeyFromSlice);
    showerItr->GetSet()->AdoptSortKeyFunc(showerKf);
    showerKf = 0;
  }

  CandTrackHandleItr * trackItr = 0;
  if (tracklist && tracklist->GetNDaughters()>0) {
    trackItr = new CandTrackHandleItr(tracklist->GetDaughterIterator());
    CandTrackHandleKeyFunc *trackKf = trackItr->CreateKeyFunc();
    trackKf->SetFun(CandTrackHandle::KeyFromSlice);
    trackItr->GetSet()->AdoptSortKeyFunc(trackKf);
    trackKf = 0;
  }

  CandSubShowerSRHandleItr * subshowerItr = 0;
  if (subshowerlist && subshowerlist->GetNDaughters()>0) {
    subshowerItr = new CandSubShowerSRHandleItr(subshowerlist->
                                                GetDaughterIterator());
    CandSubShowerSRHandleKeyFunc *subshowerKf = subshowerItr->CreateKeyFunc();
    subshowerKf->SetFun(CandSubShowerSRHandle::KeyFromSlice);
    subshowerItr->GetSet()->AdoptSortKeyFunc(subshowerKf);
    subshowerKf = 0;
  }

  Int_t nslice=0;

  // loop over slices
  CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
  while (CandSliceHandle *slice =
         dynamic_cast<CandSliceHandle*>(sliceItr())) { //loop over slices
  
    ++nslice;
    MSG("EventSS", Msg::kDebug)
      << " ****** Slice " << nslice << " *******"  
      << endl;   

    //TObjArray newevent;
    TObjArray recolist;
    TObjArray eventlist;

    //construct list of all showers and tracks in this slice
    FillRecoList(slice,showerItr,trackItr,recolist);

    //first try to remove tracks buried in showers    
    std::vector<Bool_t> removeMe(recolist.GetLast()+1);
    for (Int_t i=0;i<=recolist.GetLast();i++) {
      removeMe[i] = false;
      if(recolist.At(i)->InheritsFrom("CandTrackHandle")) {
        CandTrackHandle *track = 
          dynamic_cast<CandTrackHandle*>(recolist.At(i));
        Bool_t removeTrack = false;
        for (Int_t j=0;j<=recolist.GetLast();j++) {         
          if (recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
            CandShowerSRHandle * shower = 
              dynamic_cast<CandShowerSRHandle*>(recolist.At(j));
            if(shower->BuriedTrack(track,pMinTrackIsolation)) {
              removeTrack = true;
              break;
            }
          }
        }
        if(removeTrack){
          MSG("EventSS",Msg::kDebug) 
            << " REMOVING TRACK" << endl;
          removeMe[i] = true;
        }
      }
    }
    for (Int_t i=0;i<=recolist.GetLast();i++) {
      if(removeMe[i]) recolist.RemoveAt(i);
    }
    recolist.Compress();

    MSG("EventSS",Msg::kDebug) << "Number of reco objects in slice: " 
                               << recolist.GetLast()+1 << endl;
    //now make matrix to hold match results between all objects.
    //this also includes buried tracks that have not been removed
    //NOTE: matching requires both views match;
    //      burial just counts up matching planes in either view;
    //  One view being significantly buried may indicate a bad track, which
    //  should be associated with the shower, not used to make a new event
//gmi Change by G. Irwin 28Feb2009 -
//gmi For GCC 4.4, change old statement below to model (following that) of
//gmi http://www.codepedia.com/1/Cpp2dVector (see "One-line construction")
//gmi Original line (doesn't compile in GCC 4.4):
//gmi std::vector<std::vector<Int_t> > objectMatch(recolist.GetLast()+1,0);
//gmi Model for replacement line:
//gmi std::vector<std::vector<int> > vec(maxx,std::vector<int>(maxy,0));
//gmi Replacement line itself:
    std::vector<std::vector<Int_t> > objectMatch(recolist.GetLast()+1,
                std::vector<Int_t>(recolist.GetLast()+1,0));
    //initialize:
    for(int i=0;i<=recolist.GetLast();i++){
      for(int j=0;j<=recolist.GetLast();j++) {
//gmi Replace this line:  "objectMatch[i].push_back(0);"  by:
        objectMatch[i][j] = 0;  //gmi May be redundant with constructor init
      }
    }

    //use this vector to keep track of objects already assigned
    std::vector<Bool_t> objectUsed(recolist.GetLast()+1,0);
    
    //loop over all objects
    for(Int_t i=0;i<=recolist.GetLast();i++) { //loop over i
      objectUsed[i] = false;  //initialization
      
      CandShowerSRHandle *shower=0;
      CandTrackHandle *track=0;
      if (recolist.At(i)->InheritsFrom("CandShowerSRHandle"))
        shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
      if (recolist.At(i)->InheritsFrom("CandTrackHandle"))
        track = dynamic_cast<CandTrackHandle*>(recolist.At(i));
      
      //loop over all later objects only (since matching matrix is symmetric)
      for(Int_t j=i+1;j<=recolist.GetLast();j++) { //loop over j

        MSG("EventSS",Msg::kDebug) << "Checking objectMatch between " 
                                   << i << " and " << j << endl;

        CandShowerSRHandle *shower2=0;
        CandTrackHandle *track2=0;
        if (recolist.At(j)->InheritsFrom("CandShowerSRHandle"))
          shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(j));
        if (recolist.At(j)->InheritsFrom("CandTrackHandle"))
          track2 = dynamic_cast<CandTrackHandle*>(recolist.At(j));
        
        if(shower) {
          if(shower2) {
            if(!shower->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef) && 
               !shower2->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
              if(shower->BelongsWithShower(shower2,ac,vldcptr,pShwShwDtpos2,
                                           pShwShwDz,pShwShwDt)) objectMatch[i][j] = 1;
            }
            objectMatch[j][i] = objectMatch[i][j];
          }
          else if(track2) {
            if(shower->BelongsWithTrack(track2,ac,vldcptr,pShwTrkDtpos2,
                                        pShwTrkDz,pShwTrkDt)) objectMatch[i][j] = 1;

            if(objectMatch[i][j]!=1) {
              //check track isn't buried:
              Int_t *stats = new Int_t[8];
              shower->BuriedTrack(track2,1000,stats); //1000 forces all calculations
              if(stats[3]>0) { //at least one buried track planes
                Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
                MSG("EventSS",Msg::kDebug) 
                  << "BuriedTrack Stats: " << "\n"
                  << "Number of track planes = " << nplanes << "\n"
                  << "nSharedPlanes = " << stats[0] << "\n" 
                  << "nMissingIsoPlanes = " << stats[1] << "\n" 
                  << "nMissingSharedPlanes = " << stats[2] << "\n" 
                  << "nBuriedPlanes = " << stats[3] << "\n" 
                  << "nPhysBuriedPlanes = " << stats[4] << "\n" 
                  << "nTrkLike = " << stats[5] << "\n" 
                  << "nTrkLikeTopol = " << stats[6] << "\n" 
                  << "nXTalk = " << stats[7] << endl;
                nplanes -= (stats[1] + stats[2]); //subtract missing track planes
                //if more than half of the track planes are in the shower core: match
                if(Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac)
                  objectMatch[i][j] = -1;
                //else if more than half the track planes are anywhere in 
                //the shower or are x-talk like and are not "trklike": match
                else if(Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
                        pHardBuriedFrac) objectMatch[i][j] = -1;
                //check whether track passes through shower at high angle in either view
                //expect significant number of buried planes
                else if(stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
                        ( TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
                          TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) {
                  objectMatch[i][j] = -1;
                }
              }
              delete [] stats;
            }
            objectMatch[j][i] = objectMatch[i][j];
          }
        }
        
        if(track) {
          if(shower2) {
            if(shower2->BelongsWithTrack(track,ac,vldcptr,pShwTrkDtpos2,
                                         pShwTrkDz,pShwTrkDt)) objectMatch[i][j] = 1;
            if(objectMatch[i][j]!=1) {
              //check track isn't buried:
              Int_t *stats = new Int_t[8];
              shower2->BuriedTrack(track,1000,stats); //1000 forces all calculations
              if(stats[3]>0) { //at least one buried track planes
                Int_t nplanes = track->GetEndPlane() - track->GetBegPlane() + 1;
                MSG("EventSS",Msg::kDebug) 
                  << "BuriedTrack Stats: " << "\n"
                  << "Number of track planes = " << nplanes << "\n"
                  << "nSharedPlanes = " << stats[0] << "\n" 
                  << "nMissingIsoPlanes = " << stats[1] << "\n" 
                  << "nMissingSharedPlanes = " << stats[2] << "\n" 
                  << "nBuriedPlanes = " << stats[3] << "\n" 
                  << "nPhysBuriedPlanes = " << stats[4] << "\n" 
                  << "nTrkLike = " << stats[5] << "\n" 
                  << "nTrkLikeTopol = " << stats[6] << "\n" 
                  << "nXTalk = " << stats[7] << endl;
                nplanes -= (stats[1] + stats[2]); //subtract missing track planes
                //if more than half of the track planes are in the shower core: match
                if(Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac)
                  objectMatch[i][j] = -1;
                //else if more than half the track planes are anywhere in 
                //the shower or are x-talk like and are not "trklike": match
                else if(Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
                        pHardBuriedFrac) objectMatch[i][j] = -1;
                //check whether track passes through shower at high angle in either view
                //expect significant number of buried planes
                else if(stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
                        (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
                         TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) {
                  objectMatch[i][j] = -1;
                }
              }
              delete [] stats;
            }
            objectMatch[j][i] = objectMatch[i][j];
          }
          else if(track2) {
            if(track->BelongsWithTrack(track2,ac,vldcptr,pTrkTrkDtpos2,
                                       pTrkTrkDz,pTrkTrkDt)) objectMatch[i][j] = 1;
            objectMatch[j][i] = objectMatch[i][j];
          }
        }
        MSG("EventSS",Msg::kDebug) << "Match for objects " << i << " and " 
                                   << j << " : " << objectMatch[i][j] << endl;
      }
    }

    MSG("EventSS",Msg::kDebug) << "Starting to form events" << endl;

    //start to make events:
    Int_t nRecoObjects = recolist.GetLast()+1;
    //get total charge from each object
    std::vector<Float_t> ph(nRecoObjects,0);
    for(int i=0;i<nRecoObjects;i++){
      if (recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {
        CandShowerSRHandle *shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
        ph[i] = shower->GetCharge(CalStripType::kMIP);
      }
      if (recolist.At(i)->InheritsFrom("CandTrackHandle")) {
        CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(recolist.At(i));
        //ph[i] = track->GetCharge(CalStripType::kMIP)+1e7; //force tracks to come first
        //force tracks to come first, order tracks by number of strips
        ph[i] = track->GetNStrip()+1e7; 
      }
    }

    Int_t nEvents = 0;
    Int_t nUnusedObjects = nRecoObjects;
    while(nUnusedObjects>0) { //loop until all reco objects are in events
      
      //find largest object:
      Float_t maxPH = 0;
      Int_t maxPHIndex = -1;
      for(int i=0;i<nRecoObjects;i++){
        if(ph[i]>maxPH) {
          maxPH = ph[i];
          maxPHIndex = i;
        }
      }
      
      //add largest objects first:
      TObjArray *event = new TObjArray;
      event->Add(recolist.At(maxPHIndex));     
      ph[maxPHIndex] = 0;
      objectUsed[maxPHIndex] = true;
      nUnusedObjects-=1;

      MSG("EventSS",Msg::kDebug) << "Added object " << maxPHIndex 
                                 << " to event " << nEvents << endl;
      
      if(recolist.At(maxPHIndex)->InheritsFrom("CandTrackHandle")) {
        MSG("EventSS",Msg::kDebug) << "Object " << maxPHIndex << " is a track" << endl;

        //maxPH object is a track
        CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(recolist.At(maxPHIndex));
        
        //now add objects clearly associated with largest object:
        for(int i=0;i<nRecoObjects;i++){
          if(objectMatch[maxPHIndex][i]!=0 && !objectUsed[i]) {       
            if(recolist.At(i)->InheritsFrom("CandTrackHandle")) { 
              MSG("EventSS",Msg::kDebug) << "Track match. Adding object " << i 
                                         << " to event " << nEvents << endl;
              //if this match is a track, add it to event
              event->Add(recolist.At(i));
              ph[i] = 0;
              objectUsed[i] = true;
              nUnusedObjects-=1;
            } // end if for CandTrackHandle
            else if(recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {           
              MSG("EventSS",Msg::kDebug) << "Shower match. Testing whether to add object " 
                                         << i << " to event " << nEvents << endl;
              CandShowerSRHandle *shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
              //if this match is a shower, check where shower is relative to track vertex
              //if(TMath::Abs(shower2->GetVtxZ()-track->GetVtxZ())<pShwTrkDz/2.){
              if(TMath::Abs(shower2->GetBegPlane(PlaneView::kU) - 
                            track->GetBegPlane(PlaneView::kU)) < 
                 pShwTrkDz/(PITCHINMETRES*2.) || 
                 TMath::Abs(shower2->GetBegPlane(PlaneView::kV) - 
                            track->GetBegPlane(PlaneView::kV)) < 
                 pShwTrkDz/(PITCHINMETRES*2.)){
                MSG("EventSS",Msg::kDebug) << "Shower " << i << " is close to track vertex. "
                                           << "Adding to event " << nEvents << endl;
                event->Add(recolist.At(i));
                ph[i] = 0;
                objectUsed[i] = true;
                nUnusedObjects-=1;
                //look for objects associated with this vertex shower
                for(int j=0;j<nRecoObjects;j++){
                  if(j!=maxPHIndex && objectMatch[i][j]!=0 && !objectUsed[j]) { //got a match
                    if(objectMatch[maxPHIndex][j]==0) { // j does not match maxPHIndex
                      if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
                        MSG("EventSS",Msg::kDebug) << "Shower " << i 
                                                   << " has another track match. " 
                                                   << "Test how buried track " << j
                                                   << " is in shower" << endl; 
                        if(objectMatch[i][j]==-1){
                          event->Add(recolist.At(j));
                          ph[j] = 0;
                          objectUsed[j] = true;
                          nUnusedObjects-=1;
                        }
                        else {
                          CandTrackHandle *track2 = 
                            dynamic_cast<CandTrackHandle*>(recolist.At(j));
                          //check if track2 is buried
                          Int_t *stats = new Int_t[8];
                          shower2->BuriedTrack(track2,1000,stats);
                          if(stats[3]>0){
                            Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
                            MSG("EventSS",Msg::kDebug) 
                              << "BuriedTrack Stats: " << "\n"
                              << "Number of track planes = " << nplanes << "\n"
                              << "nSharedPlanes = " << stats[0] << "\n" 
                              << "nMissingIsoPlanes = " << stats[1] << "\n" 
                              << "nMissingSharedPlanes = " << stats[2] << "\n" 
                              << "nBuriedPlanes = " << stats[3] << "\n" 
                              << "nPhysBuriedPlanes = " << stats[4] << "\n" 
                              << "nTrkLike = " << stats[5] << "\n" 
                              << "nTrkLikeTopol = " << stats[6] << "\n" 
                              << "nXTalk = " << stats[7] << endl;
                            nplanes -= (stats[1] + stats[2]);
                            if( (Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac) ||
                                (Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
                                 pHardBuriedFrac) ||
                                (stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
                                 (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
                                  TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) ) {
                              MSG("EventSS",Msg::kDebug) << "Adding track " << j 
                                                         << " to event " << nEvents << endl;
                              event->Add(recolist.At(j));
                              ph[j] = 0;
                              objectUsed[j] = true;
                              nUnusedObjects-=1;
                            }
                          }
                          delete [] stats;
                        }
                      }
                      else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
                        MSG("EventSS",Msg::kDebug) << "Shower " << i 
                                                   << " has another shower match. " 
                                                   << "Test whether shower " << j 
                                                   << " has a better association" << endl; 
                        Bool_t isAssocTrack = false; //is this shower assoc with a track
                        for(int k=0;k<nRecoObjects;k++) { 
                          //check for match
                          if(k!=maxPHIndex && k!=i && objectMatch[j][k]!=0){
                            //check for track
                            if(recolist.At(k)->InheritsFrom("CandTrackHandle"))
                              //check it does not match with shower2
                              if(objectMatch[i][k]==0) isAssocTrack = true;
                          }
                        }
                        if(!isAssocTrack) { //no track associated
                          MSG("EventSS",Msg::kDebug) << "Shower " << j 
                                                     << " has no track association. " 
                                                     << "Add to event " << nEvents << endl; 
                          event->Add(recolist.At(j));
                          ph[j] = 0;
                          objectUsed[j] = true;
                          nUnusedObjects-=1;
                        }
                      }
                    }
                    else { //j matches maxPHIndex, so add it to event
                      MSG("EventSS",Msg::kDebug) << "Shower " << j 
                                                 << " matches orginal track " << maxPHIndex 
                                                 << " Add to event " << nEvents << endl; 
                      event->Add(recolist.At(j));
                      ph[j] = 0;
                      objectUsed[j] = true;
                      nUnusedObjects-=1;
                    }
                  }
                }
              } //end if for vertex shower
              else { //if not a vertex shower, check other associations
                MSG("EventSS",Msg::kDebug) << "Testing whether shower " << i 
                                           << " has better associations. " << endl;
                Bool_t otherTrackVertexMatch = false;
                Bool_t otherTrackMatch = false;
                Bool_t otherShowerMatch = false;
                for(int j=0;j<nRecoObjects;j++){
                  if(j!=maxPHIndex && objectMatch[i][j]!=0 && 
                     !objectUsed[j] && objectMatch[maxPHIndex][j]==0) {
                    if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
                      CandTrackHandle *track2 = 
                        dynamic_cast<CandTrackHandle*>(recolist.At(j));
                      if(TMath::Abs(shower2->GetBegPlane(PlaneView::kU) - 
                                    track->GetBegPlane(PlaneView::kU)) < 
                         pShwTrkDz/(PITCHINMETRES*2.) || 
                         TMath::Abs(shower2->GetBegPlane(PlaneView::kV) -
                                    track->GetBegPlane(PlaneView::kV)) < 
                         pShwTrkDz/(PITCHINMETRES*2.) )
                        otherTrackVertexMatch = true;
                      else {
                        //first test whether track2 is buried
                        //if so then can still associate shower to track
                        if(objectMatch[i][j]!=-1) {
                          //object match is not flagged as buried
                          //perform check:
                          Bool_t track2Buried = false;
                          Int_t *stats2 = new Int_t[8];
                          shower2->BuriedTrack(track2,1000,stats2);                       
                          if(stats2[3]>0){
                            Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
                            MSG("EventSS",Msg::kDebug) 
                              << "BuriedTrack Stats: " << "\n"
                              << "Number of track planes = " << nplanes << "\n"
                              << "nSharedPlanes = " << stats2[0] << "\n" 
                              << "nMissingIsoPlanes = " << stats2[1] << "\n" 
                              << "nMissingSharedPlanes = " << stats2[2] << "\n" 
                              << "nBuriedPlanes = " << stats2[3] << "\n" 
                              << "nPhysBuriedPlanes = " << stats2[4] << "\n" 
                              << "nTrkLike = " << stats2[5] << "\n" 
                              << "nTrkLikeTopol = " << stats2[6] << "\n" 
                              << "nXTalk = " << stats2[7] << endl;
                            nplanes -= (stats2[1] + stats2[2]);
                            if( (Float_t(stats2[4])/Float_t(nplanes) > pHardBuriedFrac) ||
                                (Float_t(stats2[3]-stats2[5]-stats2[6])/Float_t(nplanes) > 
                                 pHardBuriedFrac) ||
                                (stats2[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
                                 (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
                                  TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) ) {
                               track2Buried = true;
                              objectMatch[i][j]=-1;
                              objectMatch[j][i]=-1;
                            }
                          }
                          if(!track2Buried){
                            //test whether track2 matches better than track
                            Int_t *stats = new Int_t[8];
                            shower2->BuriedTrack(track,1000,stats);
                            if( ( stats2[3] - stats2[5] - stats2[6] ) >=
                                ( stats[3]  - stats[5]  - stats[6]  ) && 
                                ( stats2[4] >= stats[4] ) ) {
                              otherTrackMatch = true;
                            }
                            delete [] stats;
                          }
                          delete [] stats2;
                        }
                      }
                    }
                    else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
                      Int_t *stats = new Int_t[8];
                      shower2->BuriedTrack(track,1000,stats);
                      if(stats[4]==0) { //how convincing is first track/shower match?
                        otherShowerMatch = true;
                      }
                      delete [] stats;
                    }
                  }
                }
                
                //if there are no other viable matches:
                if(otherTrackVertexMatch || otherShowerMatch || 
                   otherTrackMatch) {
                  MSG("EventSS",Msg::kDebug) << "Shower " << i << " has better association: "
                                             << " otherTrackVertexMatch = " 
                                             << otherTrackVertexMatch 
                                             << " otherShowerMatch = " 
                                             << otherShowerMatch 
                                             << " otherTrackMatch = "
                                             << otherTrackMatch << endl;
                  continue;
                }
                //then add this shower to this track
                MSG("EventSS",Msg::kDebug) << "Shower " << i << " has no better association. "
                                           << "Add to event " << nEvents << endl;
                event->Add(recolist.At(i));
                ph[i] = 0;
                objectUsed[i] = true;
                nUnusedObjects-=1;                
              } //end if for non-vertex shower
            } //end if for CandShowerSRHandle
          } //end if for valid association
        } //end for loop over all objects
      } //end if maxPH object is a track
      
      else if(recolist.At(maxPHIndex)->InheritsFrom("CandShowerSRHandle")) {
        MSG("EventSS",Msg::kDebug) << "Object " << maxPHIndex << " is a shower" << endl;
        MSG("EventSS",Msg::kDebug) << "Check whether shower " << maxPHIndex 
                                   << " has other associations" << endl;
        //maxPH object is a shower
        //CandShowerSRHandle *shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(maxPHIndex));
        //now add objects clearly associated with largest object:
        for(int i=0;i<nRecoObjects;i++){
          if(objectMatch[maxPHIndex][i]!=0 && !objectUsed[i]) {   
            if(recolist.At(i)->InheritsFrom("CandTrackHandle")) {
              MSG("EventSS",Msg::kWarning) << "Arrgh - Logic flaw in eventSR. "
                                           << "Seeing tracks in event building, "
                                           << "when all tracks should already be "
                                           << "assigned!" << endl;
            } //end if CandTrackHandle
            else if(recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {
              //check for two more levels of associated showers
              //CandShowerSRHandle *shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
              MSG("EventSS",Msg::kDebug) << "Adding shower " << i << " to event " 
                                         << nEvents << endl;
              MSG("EventSS",Msg::kDebug) << "Check whether shower " << i
                                         << " has other associations" << endl;
              event->Add(recolist.At(i));
              ph[i] = 0;
              objectUsed[i] = true;
              nUnusedObjects-=1;
              for(int j=0;j<nRecoObjects;j++){
                if(j!=maxPHIndex && objectMatch[i][j]!=0 && 
                   !objectUsed[j] && objectMatch[maxPHIndex][j]==0) {
                  if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
                    MSG("EventSS",Msg::kWarning) << "Arrgh again - Logic flaw in eventSR. "
                                                 << "Seeing tracks in event building, "
                                                 << "when all tracks should already be "
                                                 << "assigned!" << endl;
                  }
                  else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
                    MSG("EventSS",Msg::kDebug) << "Adding shower " << j << " to event " 
                                               << nEvents << endl;
                    event->Add(recolist.At(j));
                    ph[j] = 0;
                    objectUsed[j] = true;
                    nUnusedObjects-=1;
                  }
                }
              }
            } //end if CandShowerSRHandle
          } //end if valid association
        } //end loop over all objects
      } //end if maxPH object is a shower
        //add event to eventlist
      eventlist.Add(event);
      nEvents+=1;
    }
    
    // remove all empty events from event list

    for (Int_t i=0; i<=eventlist.GetLast(); i++) {
      TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(i));
      if (objectlist->GetLast()<0) eventlist.Remove(objectlist);
      else if(objectlist->GetLast()==0) { //consider single shower/track events
        if(objectlist->At(0)->InheritsFrom("CandShowerSRHandle")) {
          MSG("EventSS",Msg::kDebug) << "Have a single candshower event " << endl;
          CandShowerSRHandle *showerSR = 
            dynamic_cast<CandShowerSRHandle*> (objectlist->At(0));
          //if shower seems unphysical, remove it
          if(showerSR->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
            eventlist.Remove(objectlist);
            MSG("EventSS",Msg::kDebug) << "REMOVING SHOWER-ONLY EVENT" << endl;
            continue;
          }
        }
        else if(objectlist->At(0)->InheritsFrom("CandTrackHandle")) {     
          CandTrackHandle *track = dynamic_cast <CandTrackHandle*> (objectlist->At(0));
          if(track->IsUnphysical(pTrkGapFracCut,pTrkAsymCut,pTrkXTalkFracCut,pTrkXTalkDef)) {
            //if track seems unphysical, remove it
           
            eventlist.Remove(objectlist);
            MSG("EventSS",Msg::kDebug) << "REMOVING TRACK-ONLY EVENT" << endl;
            continue;
          }
        }
      }
      else if(objectlist->GetLast()==1) { //multiple shower/track events
        for(int ind_loop = 0;ind_loop<2;ind_loop++){
          Int_t ind_loop1 = ind_loop;
          Int_t ind_loop2 = ind_loop+1; if(ind_loop2==2) ind_loop2 = 0;
          if(objectlist->At(ind_loop1)->InheritsFrom("CandShowerSRHandle")){
            CandShowerSRHandle *showerSR = 
              dynamic_cast<CandShowerSRHandle*> (objectlist->At(ind_loop1));
            if(showerSR->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
              if(objectlist->At(ind_loop2)->InheritsFrom("CandTrackHandle")){
                CandTrackHandle *track = dynamic_cast <CandTrackHandle*> (objectlist->At(ind_loop2));
                if(track->IsUnphysical(pTrkGapFracCut,pTrkAsymCut,pTrkXTalkFracCut,pTrkXTalkDef)) {
                  eventlist.Remove(objectlist);
                  MSG("EventSS",Msg::kDebug) << "REMOVING UNPHYSICAL SHOWER+TRACK EVENT" << endl;
                  continue;
                }
              }
            }
          }
        }
      }
      if(useChopper==1) {
        std::vector<CandHandle> evtVec;
        for(int i=0;i<=objectlist->GetLast();i++){
          CandHandle *event = dynamic_cast<CandHandle*>(objectlist->At(i));
          evtVec.push_back(*event);       
        }
        ChopHelper chopper(cx.GetMom());
        ChopHelp *chop = chopper.GetChopHelp(evtVec);           
        if(chop->nchop==1) continue;
        else {
          chop->Print();
          //do something with chop info
          //for(int i=0;i<ncandidates;i++){
          //}
        }
      }
    }
    eventlist.Compress();


    // now make CandEvents for each event we have found
    Bool_t buildEvent=false;
    for (Int_t i=0; i<=eventlist.GetLast(); i++) {
      TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(i));
      MSG("EventSS",Msg::kDebug)  << "making event of " 
                                    << objectlist->GetLast()+1 
                                    << " objects \n";
      cxx.SetDataIn(objectlist);
      buildEvent=true;
      CandEventHandle eventhandle = CandEvent::MakeCandidate(ah,cxx);
      MSG("EventSS",Msg::kDebug) << " # of strips:  " 
                                   << eventhandle.GetNDaughters() << "\n";
      eventhandle.SetCandSlice(slice);
      ch.AddDaughterLink(eventhandle);
      delete eventlist.At(i);
    }
    if (!buildEvent) BuildEventFromUnassoc(ac,ch,cx,slice);
  } // end loop over slice list

  delete trackItr;
  delete showerItr;
  delete subshowerItr;

  MSG("EventSS",Msg::kDebug) << "starting unassociated hits assignment \n";

  // find unassociated hits, and place them in tobjarray
  sliceItr.Reset();
  TObjArray unassociated;
  CandEventListHandle &eventlist = dynamic_cast<CandEventListHandle &>(ch);
  CandEventHandleItr * eventItr2 = 0;
  if (eventlist.GetNDaughters()>0) {    
    eventItr2 = new CandEventHandleItr(eventlist.GetDaughterIterator());
  }
  FindUnassociated(ac, sliceItr,eventItr2,unassociated);
  delete eventItr2;
  SetPrimaryShowers(ch,ac);
  
  // now loop over unnassociated hits, and check for adjacency with all
  // showers, and track vertices (in the case that no shower exists).  
  // If a match is found in the former case, add strip to shower. 
  // In the latter case, create a new shower.
 
  // this map contains the separation between each event/unassociated hit
  vector<map<const CandEventHandle*, Double_t> > dist2map(unassociated.GetLast()+1);
  FillDist2Map(ac,unassociated,ch,dist2map);
  
  // initialize a vector which holds a Bool_t
  // indicating whether a hit has already been used
  vector<Bool_t> alreadyUsed(unassociated.GetLast()+1);
  for (Int_t i=0; i<=unassociated.GetLast(); i++) alreadyUsed[i] = false;  

  // loop over unnassociated hits, adding them to events.
  // After each loop the affected distance map entries are recalculated, 
  // and the  loop is repeated.  This continues until a loop adds no strips.
  Bool_t found=true;
  Int_t n_found=0;
  while (found) {
    
    //found indicates whether a strip was added to an event in this loop
    found=false;
    
    CandStripHandle *closeststrip= 0;
    CandEventHandle *closestevent= 0;
    Double_t closestdist2 = 9999999.;
    Int_t closesti=-1;
    
    // loop over unassoc strips
    for (Int_t i=0; i<=unassociated.GetLast(); i++) {
      if(alreadyUsed[i]) continue;

      CandStripHandle *strip =
        dynamic_cast<CandStripHandle*>(unassociated.At(i));
      
      CandEventHandle *bestevent = 0;
      Double_t bestdist2 = 9999999.;
      // loop over events, and find event which is closest to strip
      // and in time
      TIter eventItr(ch.GetDaughterIterator());
      while (CandEventHandle *event =
             dynamic_cast<CandEventHandle*>(eventItr())) {
        Double_t dtime = strip->GetBegTime()-event->GetVtxT();
        if ( dtime>pHitAssocDTime0 && dtime<pHitAssocDTime1 &&
             (!bestevent || dist2map[i][event]<bestdist2) ) {
          bestevent = event;
          bestdist2 = dist2map[i][event];
        }
      }
      if (bestevent && (!closeststrip || bestdist2<closestdist2)) {
        closeststrip = strip;
        closesti=i;
        closestdist2 = bestdist2;
        closestevent = bestevent;
      }
    }
    MSG("EventSS",Msg::kVerbose)<< "closest strip " << closestdist2
                                << "/"<<pHitAssocMaxDist2<< endl;

    // if distance to closest event less than maximum allowable add to event   
    if(closeststrip && closestevent && closestdist2<pHitAssocMaxDist2){
      MSG("EventSS",Msg::kVerbose) << "closest strip "
                                   << closeststrip->GetPlane() << " "
                                   << closeststrip->GetTPos()
                                   << " " << closestdist2 << "\n";       
      alreadyUsed[closesti] = true;
      found = true;
      ++n_found;
      
      // now either add this strip to existing primary shower, or
      // create a new shower at the track vertex if one does not exist
      // Get PrimaryShower CandHandle from fShowerList for modification.
      CandShowerHandle *shower = 0;
      if (CandShowerHandle *primsh = closestevent->GetPrimaryShower()) {
        for (Int_t ishower=0; 
             !shower && ishower<closestevent->GetLastShower()+1; 
             ishower++) {
          const CandShowerHandle *target = closestevent->GetShower(ishower);
          if (primsh->IsCloneOf(*target))     // Tests clone or Cand ==
            shower = closestevent->GetShowerWritable(ishower);
        }
      }

      Float_t minShwPlane=0;
      Float_t maxShwPlane=0;
      if(shower){
        minShwPlane=min(shower->GetVtxPlane(),shower->GetEndPlane());
        maxShwPlane=max(shower->GetVtxPlane(),shower->GetEndPlane());
      }
      Float_t vertexsep = 0;
      CandTrackHandle *track = closestevent->GetPrimaryTrack();
      if (shower && track) {
        vertexsep = abs(shower->GetBegPlane() - track->GetBegPlane());
        vertexsep *= PITCHINMETRES;
        MSG("EventSS",Msg::kVerbose) << "Found track and primary shower. \n" 
                                     << "Vertex separation (limit) = " 
                                     << vertexsep << "(" 
                                     << pShwTrkDz << ")" << " \n";
      }
      // if event has shower and either no track, or shower is close
      // to track vertex, add this hit to the shower which was
      // compared in position to this strip earler.
      if (shower && (!track || vertexsep<pShwTrkDz)) {
        // NOTE: CandShowerHandle *shower, used to modify CandShower, is owned
        // by fShowerList.
        // add new strip to this shower, and the appropriate cluster
        MSG("EventSS",Msg::kVerbose) << "Shower min/max plane (limit): "
                                     << minShwPlane << "/" << maxShwPlane
                                     << " (" << pMaxNewShwLen << ")" << "\n";
        if(closeststrip->GetPlane()>minShwPlane-pMaxNewShwLen &&
           closeststrip->GetPlane()<maxShwPlane+pMaxNewShwLen ){
          MSG("EventSS",Msg::kVerbose) << "adding strip to shower \n"; 
          AddStripToEvent(closestevent,shower,subshowerlist,closeststrip);
        }
      }
      // no shower consistent with track vertex, and this strip
      // is close to vertex, so start new shower.    
      else if (track) CreatePrimaryShower(ac,cxx,closestevent,showerlist,
                                          subshowerlist,closeststrip);      
      
      // need to recalculate distances bewteen remaining
      // unassoc hits and this modified event, if the added strip
      // results in a smaller separation, replace value in dist2map 
      ReFillDist2Map(ac,closeststrip,closestevent,unassociated,dist2map);      
    } // if adding strip to event in form of new shower
  } // if found (if an unassoc hit was matched to event in last loop over them

  // because we have added strips to previously constructed showers, their
  // runalg method should be run on them again.  In the future, an
  // optimization would be to do this only on showers which have been
  // touched.
  MSG("EventSS",Msg::kVerbose) << "added " << n_found << " of "
                               << unassociated.GetLast()+1 
                               << " unassociated hits " << endl;
  
  cxx.SetDataIn(showerlist);  // Pass along persistable CandShowerList
  ReConstructShowers(ac,ch,cxx);

  //Final step: use the chopper to combine events across slices:
  if(useChopper==1) {
    std::vector<CandHandle> evtVec;
    CandEventHandleItr evitr(ch.GetDaughterIterator());
    while(CandEventHandle *event = dynamic_cast<CandEventHandle*>(evitr())){
      evtVec.push_back(*event);
    }
    ChopHelper chopper(cx.GetMom());
    ChopHelp *chop = chopper.GetChopHelp(evtVec);
    chop->Print();
  }

  // For each event, set primary shower based on 
  // distance from vertex and energy
  SetPrimaryShowers(ch,ac);
}

//____________________________________________________________________
void AlgEventSSList::AddStripToEvent(CandEventHandle * closestevent,
                                     CandShowerHandle * shower, 
                                     CandSubShowerSRListHandle * subshowerlist,
                                     CandStripHandle * closeststrip)
{
  MSG("EventSS",Msg::kVerbose) << "In AddStripToEvent:" << endl;
  closestevent->AddDaughterLink(*closeststrip);
  shower->AddDaughterLink(*closeststrip);
  if(subshowerlist && shower->InheritsFrom("CandShowerSRHandle")) {
    Calibrator& calibrator=Calibrator::Instance();
    Int_t foundsubshower=0;
    Int_t closestSubShower=-1;
    Double_t BestStripDistance = 99999;
    Double_t StripDistance = 99999;
    CandShowerSRHandle *showerSR = 
      dynamic_cast<CandShowerSRHandle*> (shower);
    MSG("EventSS",Msg::kVerbose) << "Shower UID = " 
                                 << shower->GetUidInt() << endl;

    for (Int_t isubshower=0;
         isubshower<showerSR->GetLastSubShower()+1;
         isubshower++) {
      CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
        (showerSR->GetSubShower(isubshower));
      //add strip to cluster of the same view
      if (ssl->GetPlaneView()==closeststrip->GetPlaneView()) {
        //decide which is most appropriate subshower to add strip to
        //based on distance to/spread of subshower:
        Double_t tposVtx = 0;
        if(ssl->GetPlaneView()==PlaneView::kU || 
           ssl->GetPlaneView()==PlaneView::kX) 
          tposVtx = ssl->GetVtxU();
        else if(ssl->GetPlaneView()==PlaneView::kV || 
                ssl->GetPlaneView()==PlaneView::kY)
          tposVtx = ssl->GetVtxV();
        StripDistance = 99999;
        if(ssl->GetAvgDev()>0) {
          StripDistance = 
            (TMath::Abs(closeststrip->GetTPos() - 
                        (tposVtx + ssl->GetSlope() * 
                         (closeststrip->GetZPos() - 
                          ssl->GetVtxZ()))) * 
             TMath::Cos(TMath::ATan(ssl->GetSlope()))) * 
            calibrator.GetMIP(closeststrip->
                              GetCharge(CalDigitType::kSigCorr)) / 
            ssl->GetAvgDev();
        }
        if(StripDistance>9999){ //bad slope
          StripDistance = 9999 + 
            TMath::Sqrt(TMath::Power(TMath::Abs(closeststrip->GetZPos() - 
                                                ssl->GetVtxZ()),2) +
                        TMath::Power(TMath::Abs(closeststrip->GetTPos() - 
                                                tposVtx),2));
        }
        if(StripDistance<BestStripDistance || foundsubshower==0){
          foundsubshower=1;
          BestStripDistance = StripDistance;
          closestSubShower = isubshower;
        }
      }
    }
    if(foundsubshower) {
      MSG("EventSS",Msg::kVerbose) << "adding strip to subshower "
                                   << closeststrip->GetPlane() 
                                   << endl;
      CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
        (showerSR->GetSubShower(closestSubShower));
      //extract mip -> gev conversion:
      Double_t GeVperMIP = 0;
      if(ssl->GetEnergy()!=0) GeVperMIP = ssl->GetEnergy()/
        Double_t(calibrator.GetMIP(ssl->GetCharge(CalStripType::kSigCorr)));
      TIter csshItr(subshowerlist->GetDaughterIterator());
      while (CandSubShowerSRHandle *cssh = 
             dynamic_cast<CandSubShowerSRHandle*>(csshItr())) { 
        if(cssh->IsEquivalent(ssl) && ssl->GetUidInt()==cssh->GetUidInt()) {
          MSG("EventSS",Msg::kVerbose) << "SubShower (old,new) UID " 
                                       << ssl->GetUidInt()
                                       << " " << cssh->GetUidInt() << endl; 
          showerSR->RemoveSubShower(ssl);
          CandSubShowerSRHandle *newss = cssh->DupHandle();
          newss->SetCandSlice(cssh->GetCandSlice());
          newss->AddDaughterLink(*closeststrip);
          newss->SetEnergy(GeVperMIP*Double_t(calibrator.GetMIP(closeststrip->
                 GetCharge(CalDigitType::kSigCorr))) + newss->GetEnergy());
          subshowerlist->RemoveDaughter(cssh);
          subshowerlist->AddDaughterLink(*newss);
          showerSR->AddSubShower(newss);
          delete newss;
          break;
        }
      }
    }
    else {
      MSG("EventSS",Msg::kWarning) << "Strip added to shower but not to subshower!"
                                   << " StripDistance = " << StripDistance
                                   << endl;
      MSG("EventSS",Msg::kWarning) << "Number of subshowers in shower = " 
                                   << showerSR->GetLastSubShower()+1
                                   << endl;
      for(int i=0;i<showerSR->GetLastSubShower()+1;i++){
        CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
          (showerSR->GetSubShower(i));
        MSG("EventSS",Msg::kWarning) << "SS " << i << " UID = " 
                                     << ssl->GetUidInt() << endl;
      }
    }
  }
}

//__________________________________________________________

void AlgEventSSList::CreatePrimaryShower(AlgConfig & ac, CandContext & cxx, 
                                         CandEventHandle * closestevent,
                                         CandShowerListHandle * showerlist,
                                         CandSubShowerSRListHandle * subshowerlist,
                                         CandStripHandle * closeststrip)
{
  // need to create a shower -- include strips from first plane in each
  // view plus this previously unassociated hit to generate a 3D object
  MSG("EventSS",Msg::kVerbose) << "In CreatePrimaryShower" << endl;
  // Get singleton instance of AlgFactory
  const char *pEventAlgConfig = 0;
  ac.Get("EventAlgConfig",pEventAlgConfig);

  //Double_t shwshwdz = ac.GetDouble("ShwShwDz");
  //Double_t minShwEFract = ac.GetDouble("MinShwEFract");
  Double_t pMinShwStripPE = ac.GetDouble("MinShwStripPE");
  Double_t pMaxNewShwLen = ac.GetDouble("MaxNewShwLen");

  AlgFactory &af = AlgFactory::GetInstance();
  
  CandTrackHandle *track = closestevent->GetPrimaryTrack();
  if (!showerlist) {
    MSG("EventSS",Msg::kError) <<"CandShowerListHandle missing\n" << endl;
    return;
  }
  
  closestevent->AddDaughterLink(*closeststrip);
  
  // create arrays of strips for each planeview to create
  // shower from
  
  MSG("EventSS",Msg::kVerbose) <<  "creating shower \n";
  TObjArray ustriparray;
  TObjArray vstriparray;
  Int_t uplane=track->GetBegPlane(PlaneView::kU);
  Int_t vplane=track->GetBegPlane(PlaneView::kV);
  // seed this shower with strips in most upstream track hits 
  // in each view, along with hits above a charge threshold and
  // within pMaxNewShwLen of  the track start
  TIter stripItr(track->GetDaughterIterator());
  while (CandStripHandle *tstrip =
         dynamic_cast<CandStripHandle*>(stripItr())) {
    if ((tstrip->GetPlane() == uplane) || 
        (tstrip->GetPlaneView()==PlaneView::kU &&
            tstrip->GetCharge()>pMinShwStripPE && 
            abs(tstrip->GetPlane()-uplane)<pMaxNewShwLen)) {
      ustriparray.Add(tstrip);
    }
    if ((tstrip->GetPlane() == vplane) || 
        (tstrip->GetPlaneView()==PlaneView::kV &&
            tstrip->GetCharge()>pMinShwStripPE && 
            abs(tstrip->GetPlane()-vplane)<pMaxNewShwLen)) {
      vstriparray.Add(tstrip);
    }
  }
  // add unassoc strip to appropriate array
  switch (closeststrip->GetPlaneView()) {
  case PlaneView::kU:
    ustriparray.Add(closeststrip);
    break;
  case PlaneView::kV:
    vstriparray.Add(closeststrip);
    break;
  default:
    break;
  }
  // build shower from U/V clusters or subshowers
  if(ustriparray.GetLast()>=0 && vstriparray.GetLast()>=0){
    if(showerlist->InheritsFrom("CandShowerSRListHandle") && subshowerlist){
      MSG("EventSS",Msg::kVerbose) << "U striparray has " 
                                   << ustriparray.GetLast()+1 << " entries" << endl;
      MSG("EventSS",Msg::kVerbose) << "V striparray has " 
                                   << vstriparray.GetLast()+1 << " entries" << endl;  

      // now construct CandSubShowers
      AlgHandle subshowerah = af.GetAlgHandle("AlgSubShowerSR","default");      
      cxx.SetDataIn(&ustriparray);
      CandSubShowerSRHandle usubshowerhandle =
        CandSubShowerSR::MakeCandidate(subshowerah,cxx);
      usubshowerhandle.SetCandSlice(closestevent->GetCandSlice());
      usubshowerhandle.SetClusterID(ClusterType::kUnknown);
      subshowerlist->AddDaughterLink(usubshowerhandle);
      MSG("EventSS",Msg::kVerbose) << "U subshower view: " 
                                   << usubshowerhandle.GetPlaneView() <<endl;
      
      cxx.SetDataIn(&vstriparray);
      CandSubShowerSRHandle vsubshowerhandle =
        CandSubShowerSR::MakeCandidate(subshowerah,cxx);
      vsubshowerhandle.SetCandSlice(closestevent->GetCandSlice());
      vsubshowerhandle.SetClusterID(ClusterType::kUnknown);
      subshowerlist->AddDaughterLink(vsubshowerhandle);
      MSG("EventSS",Msg::kVerbose) << "V subshower view: " 
                                   << vsubshowerhandle.GetPlaneView() <<endl;
      
      //add subshowers to list
      CandShowerSRListHandle* showerlistSR = 
        dynamic_cast<CandShowerSRListHandle*>(showerlist);
      TObjArray newshower;
      CandSubShowerSRHandle * usubshower = dynamic_cast<CandSubShowerSRHandle*>
        (subshowerlist->FindDaughter(&usubshowerhandle));
      
      CandSubShowerSRHandle * vsubshower = dynamic_cast<CandSubShowerSRHandle*>
        (subshowerlist->FindDaughter(&vsubshowerhandle));           
      newshower.Add(usubshower);
      newshower.Add(vsubshower);
      cxx.SetDataIn(&newshower);
      
      //get ShowerSS alghandle and make shower from subshowers
      AlgHandle showerahSS = af.GetAlgHandle("AlgShowerSS","default");
      CandShowerSRHandle showerhandleSR = CandShowerSR::MakeCandidate(showerahSS,cxx);
      showerhandleSR.SetCandSlice(closestevent->GetCandSlice());
      showerlistSR->AddDaughterLink(showerhandleSR);
      CandHandle *shw = showerlistSR->FindDaughter(&showerhandleSR);
      CandShowerHandle *shwh = dynamic_cast<CandShowerHandle*>(shw);
      closestevent->AddShower(shwh);
      // Change: cbs 9Aug06
      // NEED TO FORCE CREATED SHOWER TO BE PRIMARY
      // due to a change in assignment of primary shower, can now create
      // a shower here which is *not* called primary by SetPrimaryShower()
      // Causes multiple CreatePrimaryShower calls in event builder and 
      // large numbers of reco'd showers per events. 
      // SetPrimaryShowers is called at the end of RunAlg so a correct
      // assignment will be made later
      //closestevent->SetPrimaryShower(minShwEFract,shwshwdz);
      closestevent->SetPrimaryShower(shwh);
    }
  }
}

//____________________________________________________________________
void AlgEventSSList:: ReFillDist2Map(AlgConfig & ac, CandStripHandle * closeststrip,CandEventHandle * closestevent,TObjArray & unassociated,std::vector<std::map<const CandEventHandle*,Double_t> > & dist2map)
{
  Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
  Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
  Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");

  for (Int_t i=0; i<=unassociated.GetLast(); i++) {
    CandStripHandle *strip =
      dynamic_cast<CandStripHandle*>(unassociated.At(i));
    if(strip->GetPlaneView()==closeststrip->GetPlaneView()){
      Double_t dtime = strip->GetBegTime()-closestevent->GetVtxT();
      Double_t dplane =
        (Double_t)(strip->GetPlane()-closeststrip->GetPlane());
      Double_t dtpos = strip->GetTPos()-closeststrip->GetTPos();
      Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
        (pHitAssocPParm*dtpos*dtpos) +
        (pHitAssocTParm*dtime*dtime);
      if (dist2<dist2map[i][closestevent]) {
        dist2map[i][closestevent] = dist2;
      }
    }
  } // loop over unassoc hits
}

//______________________________________________________________________
void AlgEventSSList::ReConstructShowers(AlgConfig &ac, 
                                        CandHandle &ch, 
                                        CandContext &cxx)
{
  MSG("EventSS",Msg::kVerbose) << "In ReConstructShowers: " << endl;
  // Retrieve persistant modifiable CandShowerList
  CandShowerListHandle *showerlist = const_cast<CandShowerListHandle*>
    (dynamic_cast<const CandShowerListHandle*>(cxx.GetDataIn()));

  if(!showerlist->InheritsFrom("CandShowerSRListHandle")) {
    MSG("EventSS",Msg::kWarning)
      << "ListHandle is not a CandShowerSRListHandle."
      << "  ShowerList not processed in ReConstructShowers." << endl;
    return;
  }

  // Get singleton instance of AlgFactory
  AlgFactory &af = AlgFactory::GetInstance();
  
  // Get an AlgHandle to AlgShowerSR with default AlgConfig
  const char *pEventAlgConfig = 0;
  ac.Get("EventAlgConfig",pEventAlgConfig);

  AlgHandle ahss = af.GetAlgHandle("AlgShowerSS","default");
  AlgConfig acss = ahss.GetAlgConfig();
  Int_t eventCounter = 0;
  TIter evItr(ch.GetDaughterIterator());
  while (CandEventHandle *ev =
         dynamic_cast<CandEventHandle*>(evItr())) {  //event loop
    MSG("EventSS",Msg::kVerbose) << "event = " << eventCounter << endl;
    for (Int_t ishower=0; ishower<ev->GetLastShower()+1; ishower++) { //shower loop
      MSG("EventSS",Msg::kVerbose) << "shower = " << ishower << endl;
      CandShowerHandle *shower =
        dynamic_cast<CandShowerHandle*>(ev->GetShowerWritable(ishower));
      CandShowerSRHandle *showerSR = 0;
      if(shower->InheritsFrom("CandShowerSRHandle")) {
        showerSR = dynamic_cast<CandShowerSRHandle*>(shower);
      }
      MSG("EventSS",Msg::kVerbose) << "Shower UID = " << showerSR->GetUidInt() << endl;
      //assume subshower chain
      if(showerSR) {  
        TObjArray newshower;
        // fill array of subshowers for this shower
        for (Int_t isubshower=0;
             isubshower<showerSR->GetLastSubShower()+1;
             isubshower++) {
          MSG("EventSS",Msg::kVerbose) << "subshower = " << isubshower << endl;
          const CandSubShowerSRHandle *ssh = 
            showerSR->GetSubShower(isubshower);
          if (ssh) {
            // New handle on heap
            CandSubShowerSRHandle *sshd = ssh->DupHandle();
            // Add new heap-based CandHandle to TObjArray
            newshower.Add(sshd);
            MSG("EventSS",Msg::kDebug) << "Uid's (old,new) " 
                                       << ssh->GetUidInt() << " " 
                                       << sshd->GetUidInt() << endl;
          }
        }
        // re-run algorithm for this shower
        if (newshower.GetEntriesFast() > 0) {
          cxx.SetDataIn(&newshower);
          ahss.RunAlg(*showerSR,cxx);
        }
        else {
          MSG("EventSS", Msg::kError)
            << "Attempt to pass empty TObjArray newshower to AlgShowerSS"
            << " in CandContext.  AlgShowerSS::RunAlg() call ignored."
            << endl;
        }
        // Delete heap-based CandHandles in TObjArray
        newshower.Delete();
      }
      else {
        MSG("EventSS",Msg::kWarning)
          << "Handle is not a CandShowerSRHandle."
          << "  Shower not processed in ReConstructShowers."
          << endl;
        continue;
      }
      MSG("EventSS",Msg::kVerbose) << "Modified shower UID = " 
                                   << showerSR->GetUidInt() << endl; 
      // now loop through strips in this shower, and determine whether primary
      // track in this event intercepts this strip.  If so, subtract 1 mip
      // from shower energy.
      
      // if event has primary track, extrapolate through shower 
      // and subtract energy loss to recalibrate 
      CandTrackHandle *primaryTrack = ev->GetPrimaryTrack();
      if (primaryTrack) {
        shower->CalibrateEnergy(primaryTrack,acss);
      }
      
      // Update modified showers in persistent modifiable CandShowerList.
      // CandHandle could potentially supply a ReplaceDaughter() function.
      // Interface of ReplaceDaughter() would be debatable.
      if (showerlist) {
        TIter shiter(showerlist->GetDaughterIterator());
        CandShowerHandle *target;
        Bool_t found = kFALSE;
        while (!found &&
               (target = dynamic_cast<CandShowerHandle*>(shiter()))) {
          if (shower->IsCloneOf(*target)) {  // Tests clone or Cand ==
            found = kTRUE;
            
            // Replace target only if shower is a modified version of target
            if (!(shower->IsEqual(target))) {   // CandBase inequality
              shower->SetCandSlice(ev->GetCandSlice());
              if (!showerlist->RemoveDaughter(target)) {    // Failure
                MSG("EventSS",Msg::kError) 
                  << endl
                  << "Failure of ShowerList daughter removal " << endl
                  << "attempted during replacement by modified Shower."
                  << endl << "Will result in double counted Shower."
                  << endl;
              }
              showerlist->AddDaughterLink(*shower);
            }
          }
        }
        
        // Add new (unfound) shower to persistent modifiable CandShowerList.
        if (!found){
          shower->SetCandSlice(ev->GetCandSlice());
          showerlist->AddDaughterLink(*shower);
        }
      }
    }
    eventCounter++;
  }
}

//______________________________________________________________________
void AlgEventSSList::SetPrimaryShowers(CandHandle &ch,AlgConfig &ac)
{
  
  // This method loops over events in the eventlist, finding the largest
  // shower and setting that to the primary shower in the event.
  
  Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
  Double_t pMinShwEFract = ac.GetDouble("MinShwEFract");
  
  TIter evItr(ch.GetDaughterIterator());
  while (CandEventHandle *ev =
                              dynamic_cast<CandEventHandle*>(evItr())) {
    ev->SetPrimaryShower(pMinShwEFract,pShwShwDz);
  }
}

//______________________________________________________________________
void AlgEventSSList::BuildEventFromUnassoc(AlgConfig &ac,
                                           CandHandle &ch,
                                           CandContext &cx, 
                                           CandSliceHandle * slice)
{
  assert(cx.GetDataIn());
  MSG("EventSS",Msg::kDebug) <<" Building Event - no showers or tracks " 
                             << endl;
  if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
    MSG("EventSS",Msg::kDebug)<<" dataIn inherits from TObjArray " << endl;
    return;
  }

  Double_t pShwXTalkFracCut = ac.GetDouble("ShwXTalkFracCut");
  Double_t pShwXTalkDef = ac.GetDouble("ShwXTalkDef");

  const CandSliceListHandle *slicelist = 0;
  CandShowerListHandle *showerlist = 0;
  CandSubShowerSRListHandle *subshowerlist = 0;
  const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
  for (Int_t i=0; i<=cxin->GetLast(); i++) {
    TObject *tobj = cxin->At(i);
    if (tobj->InheritsFrom("CandSliceListHandle")) {
      slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
    }
    // Retrieve persistant modifiable CandShowerList
    if (tobj->InheritsFrom("CandShowerListHandle")) {
      showerlist = dynamic_cast<CandShowerListHandle*>(tobj);
    }
    if (tobj->InheritsFrom("CandSubShowerSRListHandle")) {
      subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj);
    }
  }
  
  if (!showerlist || !slicelist || !subshowerlist) {
    MSG("EventSS",Msg::kDebug)
      << "Missing either slicelist, showerlist or subshowerlist" << endl;
    return;
  }
  Int_t detectorType = slicelist->GetVldContext()->GetDetector();
  
  CandContext cxx(this,cx.GetMom());
  AlgFactory &af = AlgFactory::GetInstance();
  
  const char *pEventAlgConfig = 0;
  ac.Get("EventAlgConfig",pEventAlgConfig);
  AlgHandle ah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
  
  Double_t pHitAssocDPlane = ac.GetInt("HitAssocDPlane");
  Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
  Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
  Double_t pHitAssocMaxDist2 = ac.GetDouble("HitAssocMaxDist2");
  Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
  Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
  Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");
  Double_t pShwShwDtpos2 = ac.GetDouble("ShwShwDtpos2");
  Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
  Double_t pShwShwDt = ac.GetDouble("ShwShwDt");
  Double_t pMinShwEFract = ac.GetDouble("MinShwEFract");
  Double_t pMinUnassocStripPE = ac.GetDouble("MinUnassocStripPE");


  const CandRecord *candrec = cx.GetCandRecord();
  assert(candrec);
  const VldContext *vldcptr = candrec->GetVldContext();
  assert(vldcptr);
  VldContext vldc = *vldcptr;
  UgliGeomHandle ugh(vldc);

  // Iterate over strips , finding sets which are 'in time', and within
  // a given distance of the nearest neighbor in the set, in both views.
  // Each of these sets is then used to construct a shower, which becomes
  // the primary shower of a new event.  We continue until an iteration
  // in which no new event is created.
  
  // fill array of unique unassociated strips
  TObjArray unassociated;
  TIter stripItr(slice->GetDaughterIterator());
  while (CandStripHandle *strip =
         dynamic_cast<CandStripHandle*>(stripItr())) {

    // Ignore strips below the required minimum PE
    if (strip->GetCharge(CalDigitType::kPE)<pMinUnassocStripPE) continue;

    Bool_t found=false;
    for (Int_t i=0; i<=unassociated.GetLast() && !found; i++) {
      CandStripHandle *strip0 =
        dynamic_cast<CandStripHandle*>(unassociated.At(i));
      if (*strip == *strip0) {
        found = true;
      }
    }
    if (!found) {
      if(!(detectorType==Detector::kNear && strip->GetPlane()>120)){
        unassociated.Add(strip);
      }
    }
  }
  
  // loop until no new matches are found in the loop over unassoc hits.
  Bool_t found=false;
  while (unassociated.GetLast()>0 && !found) {
    MSG("EventSS",Msg::kDebug)<< " # unassoc " 
                              << unassociated.GetLast()+1 << endl; 
    Bool_t firstu=true;
    Bool_t firstv=true;
    TObjArray neweventu;
    TObjArray neweventv;
    Bool_t addedstrip = true;
    Float_t totcharge=0;
    while (addedstrip) {
      addedstrip=false;      
      // loop over unassoc strips.
      for (Int_t i=0; i<=unassociated.GetLast(); i++) {
        CandStripHandle *strip =
          dynamic_cast<CandStripHandle*>(unassociated.At(i));
        switch (strip->GetPlaneView()) {
        case PlaneView::kU:
          // strip is in U plane
          // if this is the first strip, add to newevent strip list
          if (firstu) {
            // if we don't have v strips yet, just add this strip to u list
            // seed the event
            if (firstv) {
              neweventu.Add(strip);
              MSG("EventSS",Msg::kDebug)<< " add  first  " << endl;
              totcharge +=strip->GetCharge();
              addedstrip = true;
              firstu=false;
            }
            // if we already have one or more v strips check that they are
            // in time, and  agree in Z position.
            else {
              Bool_t found2=false;
              for (Int_t j=0; j<=neweventv.GetLast() && !found2;
                   j++) {
                CandStripHandle *estrip =
                  dynamic_cast<CandStripHandle*>(neweventv.At(j));
                Double_t dtime =
                  strip->GetBegTime()-estrip->GetBegTime(); 
                MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
                                             << " " <<pHitAssocDTime0*1e9 
                                             << "/" << pHitAssocDTime1*1e9 
                                             << " plane diff " 
                                             <<  abs(strip->GetPlane() - 
                                                     estrip->GetPlane()) 
                                             << "/" << pHitAssocDPlane << endl;
                if (dtime>pHitAssocDTime0 &&
                    dtime<pHitAssocDTime1 &&
                    abs(strip->GetPlane()-estrip->GetPlane()) <=
                    pHitAssocDPlane) {
                  MSG("EventSS",Msg::kVerbose)<< " add  " << endl;
                  // match is found, add this strip to array
                  neweventu.Add(strip);
                  totcharge +=strip->GetCharge();
                  addedstrip = true;
                  firstu=0;
                  found2=true;
                } // if match found
              }  // for loop over V strips in new event
            } // end not first V strip
          } // end first U strip
          else {  
            // if we already have u strips, make sure that the one we
            // add is adjaccent, using distance criteria based
            // on plane, time, and transverse postion separation
            Bool_t found2=false;
            for (Int_t j=0; j<=neweventu.GetLast() && !found2; j++) {
              CandStripHandle *estrip =
                dynamic_cast<CandStripHandle*>(neweventu.At(j));
              Double_t dtime =
                strip->GetBegTime()-estrip->GetBegTime(); 
              MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
                                           << " " <<pHitAssocDTime0*1e9 
                                           << "/" << pHitAssocDTime1*1e9 
                                           << " plane diff " 
                                           << abs(strip->GetPlane() - 
                                                  estrip->GetPlane())
                                           << " dtpos " << (strip->GetTPos() - 
                                                            estrip->GetTPos())
                                           << endl;
              if (dtime>pHitAssocDTime0 &&
                  dtime<pHitAssocDTime1) {
                Double_t dplane =
                  (Double_t)(strip->GetPlane()-estrip->GetPlane());
                Double_t dtpos = strip->GetTPos()-estrip->GetTPos();
                Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
                  (pHitAssocPParm*dtpos*dtpos) +
                  (pHitAssocTParm*dtime*dtime);
                // if dist2 2 parameter meets match criteria
                // add strip to array
                if (dist2<pHitAssocMaxDist2*4) {
                  MSG("EventSS",Msg::kVerbose)<< " add " << endl;
                  neweventu.Add(strip);
                  totcharge +=strip->GetCharge();
                  addedstrip = true;
                  firstu=0;
                  found2=true;
                } // if positions match
              } // if in time
            } // end loop over event u strips
          } // end if already have u strip
          break;
          
        case PlaneView::kV:
          
          // if this is the first strip, add to newevent strip list
          if (firstv) {
            // strip is in V plane
            // if this is the first strip, add to newevent strip list
            if (firstu) {
              neweventv.Add(strip);
              totcharge +=strip->GetCharge();
              MSG("EventSS",Msg::kDebug)<< " add  first" << endl;
              addedstrip = true;
              firstv = 0;
            }
            else {
              // if we already have one or more u strips check that they are
              // in time, and require agreement in Z position.
              Bool_t found2=false;
              for (Int_t j=0; j<=neweventu.GetLast() && !found2;
                   j++) {
                CandStripHandle *estrip =
                  dynamic_cast<CandStripHandle*>(neweventu.At(j));
                Double_t dtime =
                  strip->GetBegTime()-estrip->GetBegTime();
                MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
                                             << " " << pHitAssocDTime0*1e9 
                                             << "/" << pHitAssocDTime1*1e9 
                                             << " plane diff " 
                                             <<  abs(strip->GetPlane() - 
                                                     estrip->GetPlane()) 
                                             << "/" << pHitAssocDPlane << endl;
                if (dtime>pHitAssocDTime0 &&
                    dtime<pHitAssocDTime1 &&
                    abs(strip->GetPlane()-estrip->GetPlane()) <=
                    pHitAssocDPlane) {
                  // if dist2 2 parameter meets match criteria, 
                  //add strip to array
                  MSG("EventSS",Msg::kVerbose) << " add  " << endl;
                  neweventv.Add(strip);
                  totcharge +=strip->GetCharge();
                  addedstrip = true;
                  firstv = false;
                  found2 = true;
                }
              } // end loop over u strips in event
            } // end if already have u strips
          } // end if this is the first v strip
          else {
            // if we already have v strips, make sure that the one we
            // add is adjaccent, using same distance criteria base
            // on plane, time, and transverse postion separation
            Bool_t found2=false;
            for (Int_t j=0; j<=neweventv.GetLast() && !found2; j++) {
              CandStripHandle *estrip =
                dynamic_cast<CandStripHandle*>(neweventv.At(j));
              Double_t dtime = strip->GetBegTime() -
                estrip->GetBegTime();              
              MSG("EventSS",Msg::kVerbose)<< " dtime " << dtime*1e9 << " " 
                                          <<  pHitAssocDTime0*1e9 << "/" 
                                          << pHitAssocDTime1*1e9 
                                          << " plane diff " 
                                          <<  abs(strip->GetPlane() - 
                                                  estrip->GetPlane())
                                          << " dtpos " << (strip->GetTPos() - 
                                                           estrip->GetTPos())
                                          << endl;
              if (dtime>pHitAssocDTime0 &&
                  dtime<pHitAssocDTime1) {
                Double_t dplane =
                  (Double_t)(strip->GetPlane()-estrip->GetPlane());
                Double_t dtpos = strip->GetTPos()-estrip->GetTPos();
                Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
                  (pHitAssocPParm*dtpos*dtpos) +
                  (pHitAssocTParm*dtime*dtime);
                // if dist2 2 parameter meets match criteria, 
                // add strip to array
                if (dist2<pHitAssocMaxDist2*4) {
                  MSG("EventSS",Msg::kVerbose) << " add  " << endl;
                  neweventv.Add(strip);
                  totcharge +=strip->GetCharge();
                  addedstrip = true;
                  firstv = false;
                  found2 = true;
                } // end if add strip
              } //end if in time
            }  // end loop over v strips
          }
          break;
          
        default:
          break;
        }
      }
      
      // end of loop over hits - remove new 'associations' and compress
      for (Int_t i=0; i<=neweventu.GetLast(); i++) {
        CandStripHandle *estrip =
          dynamic_cast<CandStripHandle*>(neweventu.At(i));
        unassociated.Remove(estrip);
      }
      for (Int_t i=0; i<=neweventv.GetLast(); i++) {
        CandStripHandle *estrip =
          dynamic_cast<CandStripHandle*>(neweventv.At(i));
        unassociated.Remove(estrip);
      }
      unassociated.Compress();
    }
    
    // if we have hits in both views, build an event from them
    // if pulse height>10 pes
    MSG("EventSS",Msg::kDebug)<< " considering event formation  u/v: " 
                              << !firstu  << "/" << !firstv << " charge " 
                              << totcharge << endl;
    if (!firstu && !firstv && totcharge>20.0) {      
      found=true;

      //construct CandSubShowerSRs
      cxx.SetDataIn(&neweventu);
      AlgHandle subshowerah = af.GetAlgHandle("AlgSubShowerSR","default");
      CandSubShowerSRHandle usubshowerhandle =
        CandSubShowerSR::MakeCandidate(subshowerah,cxx);
      usubshowerhandle.SetCandSlice(slice);
      usubshowerhandle.SetClusterID(ClusterType::kUnknown);
      subshowerlist->AddDaughterLink(usubshowerhandle);
      cxx.SetDataIn(&neweventv);
      CandSubShowerSRHandle vsubshowerhandle =
        CandSubShowerSR::MakeCandidate(subshowerah,cxx);
      vsubshowerhandle.SetCandSlice(slice);
      vsubshowerhandle.SetClusterID(ClusterType::kUnknown);
      subshowerlist->AddDaughterLink(vsubshowerhandle);
      CandHandle *usubshower=subshowerlist->FindDaughter(&usubshowerhandle);
      CandHandle *vsubshower=subshowerlist->FindDaughter(&vsubshowerhandle);
      
      //build shower
      TObjArray newshower;
      newshower.Add(usubshower);
      newshower.Add(vsubshower);
      cxx.SetDataIn(&newshower);
      AlgHandle showerah = af.GetAlgHandle("AlgShowerSS","default");
      CandShowerSRHandle showerhandle =
        CandShowerSR::MakeCandidate(showerah,cxx);
      showerhandle.SetCandSlice(slice);
      showerhandle.SetCandRecord(slicelist->GetCandRecord());

      Bool_t newShowerOverlapsOld=false;
      CandShowerHandleItr showerItr(showerlist->GetDaughterIterator());      
      while (CandShowerHandle *shower =
             dynamic_cast<CandShowerHandle*>(showerItr())){

        Double_t dz = showerhandle.GetVtxZ()-shower->GetVtxZ();
        Double_t du = showerhandle.GetVtxU()-shower->GetVtxU();
        Double_t dv = showerhandle.GetVtxV()-shower->GetVtxV();
        Double_t dt = showerhandle.GetVtxT()-shower->GetVtxT();

        if( ( (du*du+dv*dv<pShwShwDtpos2 &&fabs(dz)<pShwShwDz) ||
              (du*du<pShwShwDtpos2/4. && dv*dv<pShwShwDtpos2*4. && 
               fabs(dz)<pShwShwDz/2.) ||
              (du*du<pShwShwDtpos2*4 && dv*dv<pShwShwDtpos2/4. && 
               fabs(dz)<pShwShwDz/2.) ) &&  
            fabs(dt)<pShwShwDt ) {
          newShowerOverlapsOld = true;
          MSG("EventSS",Msg::kDebug) 
            << " new shower overlaps with previous - don't make new event " 
            << endl;
          break;
        }
      }
        
      if(!newShowerOverlapsOld){
        if(!showerhandle.IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
          MSG("EventSS",Msg::kDebug)<< " Building Event! " << endl;
          showerlist->AddDaughterLink(showerhandle);
          TObjArray recolist;
          recolist.Add(&showerhandle);
          cxx.SetDataIn(&recolist);
          AlgHandle eventah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
          CandEventHandle eventhandle =
            CandEvent::MakeCandidate(eventah,cxx);
          eventhandle.SetCandSlice(slice);
          eventhandle.SetPrimaryShower(pMinShwEFract,pShwShwDz);
          ch.AddDaughterLink(eventhandle);      
        }
      }
    }
  }
}

//_________________________________________________________________________________
// for a given slice, fill an array with all tracks and showers
void AlgEventSSList::FillRecoList(CandSliceHandle *slice, CandShowerHandleItr * showerItr, 
                                  CandTrackHandleItr* trackItr, TObjArray & recolist){

  if (showerItr) {
    showerItr->Reset();
    while (CandShowerHandle *shower =
           dynamic_cast<CandShowerHandle*>((*showerItr)())) {
      if (shower->GetCandSlice()) {
        if (slice->IsCloneOf(*(shower->GetCandSlice()))) {
          recolist.Add(shower);
        }
      }
      else {
        MSG("EventSS", Msg::kError)
          << "Shower doesn't contain a slice.  Not added to recolist."
          << endl;
      }
    }
  }
  
  if (trackItr) {
    trackItr->Reset();
    while (CandTrackHandle *track =
           dynamic_cast<CandTrackHandle*>((*trackItr)())) {
      if (track->GetCandSlice()) {
        if (slice->IsCloneOf(*track->GetCandSlice())) {
          if(track->GetNDaughters()==0 && 
             track->InheritsFrom("CandFitTrackHandle") &&  
             dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack()){
            track = dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack();    
            MSG("EventSS", Msg::kDebug)<< "Finder Track being used " << endl;
          }
          recolist.Add(track);
        }
      }
      else {
        MSG("EventSS", Msg::kError)
          << "Track doesn't contain a slice.  Not added to recolist."
          << endl;
      }
    }
  }
}

//_______________________________________________________________
void  AlgEventSSList::AddObjectToEvent(CandRecoHandle * reco, TObjArray *objectlist, 
                                       TObjArray * prevlist, Bool_t merge ){
  if (!merge) {
    // if this is the first match we have found for this object, 
    //add object to list
    // and point prevlist to this object list
    objectlist->Add(reco);
    MSG("EventSS",Msg::kVerbose) << "    adding to list\n";
  }
  else {                  
    MSG("EventSS",Msg::kVerbose) << "    combining lists\n";
    // this object matches more than one event - merge their member
    // objects to list prevlist, and add the new object
    for (Int_t ireco2=0; ireco2<=objectlist->GetLast();ireco2++)
      prevlist->Add(objectlist->At(ireco2));
    prevlist->Add(reco);
    objectlist->Clear();
  }
}

//__________________________________________________________________________
// fill TOjArray with strips not in any object daughter list in this snarl
void AlgEventSSList::FindUnassociated(AlgConfig& ac, CandSliceHandleItr & sliceItr, 
                                      CandEventHandleItr *eventItr, 
                                      TObjArray & unassociated){

  Double_t pMinUnassocStripPE = ac.GetDouble("MinUnassocStripPE");

  while (CandSliceHandle *slice =
         dynamic_cast<CandSliceHandle*>(sliceItr())) {
    TIter stripItr(slice->GetDaughterIterator());
    while (CandStripHandle *strip =
           dynamic_cast<CandStripHandle*>(stripItr())) {

      // Ignore strips below the required minimum PE
      if (strip->GetCharge(CalDigitType::kPE)<pMinUnassocStripPE) continue;

      Bool_t found=false;
      // check for duplicates      
      for (Int_t i=0; i<=unassociated.GetLast() && !found; i++) {
        CandStripHandle *strip0 =
          dynamic_cast<CandStripHandle*>(unassociated.At(i));
        if (*strip == *strip0) {
          found = true;
          break;
        }
      }
      // if this strip not already in unnassociated list, then
      // loop over objects in this slice, and check against all daughter strips
      if (!found) {
        if (eventItr) {
          eventItr->Reset();
          while (CandEventHandle *event =
                 dynamic_cast<CandEventHandle*>((*eventItr)())) {
            if (event->GetCandSlice()) {
              if (slice->IsCloneOf(*(event->GetCandSlice()))) { 
                if (event->FindDaughter(strip)) {
                  found = true;
                  break;
                }
              }
            }
          }
        }
      }
      // this strip is not in any event, add to unassociated list
      if (!found) {
        unassociated.Add(strip);
      }
    }
  }
}

//____________________________________________________________________

void AlgEventSSList::FillDist2Map(AlgConfig & ac, TObjArray & unassociated,CandHandle & ch,vector<map<const CandEventHandle*, Double_t> > & dist2map){
  // grab all needed algorithm parameters  
  
  Double_t pShwTrkDz = ac.GetDouble("ShwTrkDz");
  Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
  Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
  Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
  Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
  Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");

  // construct the map, first looping over strips
  for (Int_t i=0; i<=unassociated.GetLast(); i++) {
    CandStripHandle *strip =
      dynamic_cast<CandStripHandle*>(unassociated.At(i));
    TIter eventItr(ch.GetDaughterIterator());
    // inner loop is over events
    while (CandEventHandle *event =
           dynamic_cast<CandEventHandle*>(eventItr())) {
      
      // calculate separation in space/time
      Double_t dtime = (strip->GetBegTime()-event->GetVtxT());
      Double_t bestdist2=-1.;
      Double_t bestdplane=-1;
      Double_t bestdtpos=-1;
      Double_t bestdtime=-1;
      dist2map[i][event] = 999999999.;
      MSG("EventSS",Msg::kVerbose) << i << " " << strip->GetPlane() << " " 
                                   << strip->GetStrip() << " stp beg time "
                                   << strip->GetBegTime() << " evt vtx time "
                                   << event->GetVtxT()
                                   << " dtime " 
                                   << dtime*1e9 << "/" << pHitAssocDTime0*1e9 
                                   << "/" << pHitAssocDTime1*1e9 << endl;  
      // if there is a rough time match proceed
      if (dtime>pHitAssocDTime0 && dtime<pHitAssocDTime1) {

        Float_t vertexsep = 0;
        // If  event has both shower
        // and track calc separation between the two.  If this 
        // separation is greater
        // than pSHwTrkDz, we won't bother adding the strip to the
        // primary shower, and will instead consider whether
        // to create a new shower at track vertex.
        if (event->GetPrimaryShower() && event->GetPrimaryTrack())
          vertexsep = fabs(event->GetPrimaryShower()->GetVtxZ() -
                          event->GetPrimaryTrack()->GetVtxZ());

        // the event has a shower, and either no track, or track and
        // shower vertices are in agreement.  In this case, we
        // consider adding this strip to this shower, so calc separation.
        if (event->GetPrimaryShower() && 
               (!event->GetPrimaryTrack() || vertexsep<pShwTrkDz)) {

          // create iterator of shower daughter strips in same view
          // as unassoc strip
          CandStripHandleItr shwstripItr(event->GetPrimaryShower()->
                                         GetDaughterIterator());
          CandStripHandleKeyFunc *shwstripKf = shwstripItr.CreateKeyFunc();
          shwstripKf->SetFun(CandStripHandle::KeyFromView);
          shwstripItr.GetSet()->AdoptSortKeyFunc(shwstripKf);
          shwstripKf = 0;
          shwstripItr.GetSet()->Slice();
          shwstripItr.GetSet()->Slice(strip->GetPlaneView(),
                                      strip->GetPlaneView());
          
          // iterate over the daughter list, and find minimum distane
          // between daughter and unassoc strip
          while (CandStripHandle *shwstrip =
                        dynamic_cast<CandStripHandle*>(shwstripItr())) {
            Double_t dplane = (Double_t)(strip->GetPlane()-shwstrip->GetPlane());
            Double_t dtpos = strip->GetTPos()-shwstrip->GetTPos();
            Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
                             (pHitAssocPParm*dtpos*dtpos) +
                             (pHitAssocTParm*dtime*dtime);
            if (bestdist2<0. || dist2<bestdist2) {
              bestdist2 = dist2;
              bestdplane=dplane;
              bestdtpos=dtpos;
              bestdtime=dtime;
            }
          }
          MSG("EventSS",Msg::kVerbose)
            << "primary shower." << " dplane:" << bestdplane <<" dtpos:"
            << bestdtpos <<" dtime:" << bestdtime <<" dist2:" << bestdist2<< "\n";
          dist2map[i][event] = bestdist2;
        } // end if event has shower, etc.
        // else if criteria above was not satisfied, consider whether
        // strip is near track vertex, in which case we start new shower
        else if (event->GetPrimaryTrack()) {

          // calc distance to vertex
          Double_t dplane = (Double_t)(strip->GetPlane() - 
                                       event->GetPrimaryTrack()->
                                       GetBegPlane());
          Double_t dtpos=0;
          switch (strip->GetPlaneView()) {
          case PlaneView::kU:
            dtpos=strip->GetTPos()-event->GetPrimaryTrack()->
              GetU(event->GetPrimaryTrack()->GetBegPlane());
            break;
          case PlaneView::kV:
            dtpos=strip->GetTPos()-event->GetPrimaryTrack()->
              GetV(event->GetPrimaryTrack()->GetBegPlane());
            break;
          default:
            break;
          }
          Double_t dist2 = ( (pHitAssocZParm*dplane*dplane) +
                             (pHitAssocPParm*dtpos*dtpos) +
                             (pHitAssocTParm*dtime*dtime) );
          dist2map[i][event] = dist2;
          MSG("EventSS",Msg::kVerbose)
            << "primary track. begin:" << " dplane:" << dplane
            << " dtpos:" << dtpos <<" dtime:" << dtime <<" dist2:"
            << dist2 << "\n";
        } // end if primary track
      }  // end if time match OK
    }  // end loop over events
  }  // end loop over unassoc strips
}

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

---