AlgEventSRList.cxx

Go to the documentation of this file.

// $Id: AlgEventSRList.cxx,v 1.67 2008/11/05 17:26:55 rodriges Exp $
//
// AlgEventSRList.cxx
//
// AlgEventSRList is a concrete EventSRList 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/AlgEventSRList.h"
#include "RecoBase/CandEvent.h"
#include "RecoBase/CandEventHandle.h"
#include "RecoBase/CandEventList.h"
#include "RecoBase/CandEventListHandle.h"
#include "Conventions/PlaneView.h"
#include "RecoBase/CandClusterHandle.h"
#include "RecoBase/CandCluster.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 "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandClusterListHandle.h"
#include "RecoBase/CandSliceListHandle.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"
ClassImp(AlgEventSRList)

#define PITCHINMETRES 0.0594

//______________________________________________________________________
CVSID("$Id: AlgEventSRList.cxx,v 1.67 2008/11/05 17:26:55 rodriges Exp $");

//______________________________________________________________________
AlgEventSRList::AlgEventSRList()
{
}

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

//______________________________________________________________________
void AlgEventSRList::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  
  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");

// 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;
  CandClusterListHandle *clusterlist = 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("CandClusterListHandle")) {
      clusterlist = dynamic_cast<CandClusterListHandle*>(tobj);
    }
  }
  if (!slicelist) {
    MSG("EventSR",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;
  }

  CandClusterHandleItr * clusterItr = 0;
  if (clusterlist && clusterlist->GetNDaughters()>0) {
    clusterItr = new CandClusterHandleItr(clusterlist->GetDaughterIterator());
    CandClusterHandleKeyFunc *clusterKf = clusterItr->CreateKeyFunc();
    clusterKf->SetFun(CandClusterHandle::KeyFromSlice);
    clusterItr->GetSet()->AdoptSortKeyFunc(clusterKf);
    clusterKf = 0;
  }

  Int_t nslice=0;

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

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

    // remove tracks in recolist which have few/any isolated hits
    RemoveTracksinShowers(ac,&recolist);
    
    MSG("EventSR",Msg::kVerbose)
      << "recolist begin " << recolist.GetLast()+1 << "\n";
    // loop over shower+tracklist, and find associations
    for (Int_t i=0; i<=recolist.GetLast(); i++) {
      //determine whether this entry is shower or track
      CandShowerHandle *shower=0;
      CandTrackHandle *track=0;
      CandRecoHandle *reco=0;
      if (recolist.At(i)->InheritsFrom("CandShowerHandle")) {
        shower = dynamic_cast<CandShowerHandle*>(recolist.At(i));
        reco = shower;
        MSG("EventSR",Msg::kVerbose) << "  shower " << reco->GetVtxU()
                                     << " " << reco->GetVtxV() << " " 
                                     << reco->GetVtxZ() << "\n";
      }
      if (recolist.At(i)->InheritsFrom("CandTrackHandle")) {
        track = dynamic_cast<CandTrackHandle*>(recolist.At(i));
        reco = track;
        MSG("EventSR",Msg::kVerbose) << "  track  " << reco->GetVtxU()
                                     << " " << reco->GetVtxV() << " " 
                                     << reco->GetVtxZ() << "\n";
      }
      
      // now loop over events previously created
      // and see whether this object should be added
      Bool_t merge=false;
      TObjArray *prevlist=0;
      for (Int_t ievt=0; ievt<=eventlist.GetLast(); ievt++) {
        Bool_t thiseventmatches=false;
        TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(ievt));
        // loop over objects within this event
        for (Int_t ireco=0; !thiseventmatches && ireco<=objectlist->GetLast(); ireco++) {
          Bool_t match=false;
          if (objectlist->At(ireco)->InheritsFrom("CandShowerHandle")) {
            CandShowerHandle *target = dynamic_cast<CandShowerHandle*>
              (objectlist->At(ireco));
            if(shower) {
              match = target->BelongsWithShower(shower,ac,vldcptr,
                                                pShwShwDtpos2,
                                                pShwShwDz,pShwShwDt);
            }
            if(track) {
              match = target->BelongsWithTrack(track,ac,vldcptr,
                                               pShwTrkDtpos2,
                                               pShwTrkDz,pShwTrkDt);
              match = match && MatchOtherTracksInEvent(ac,vldcptr,track,objectlist);
            }
          }
          else if (objectlist->At(ireco)->InheritsFrom("CandTrackHandle")) {
            CandTrackHandle *target = dynamic_cast<CandTrackHandle*>
              (objectlist->At(ireco));
            if(shower) {
              match = target->BelongsWithShower(shower,ac,vldcptr,
                                                pShwTrkDtpos2,
                                                pShwTrkDz,pShwTrkDt);         
            }
            if(track) {
              match = target->BelongsWithTrack(track,ac,vldcptr,
                                               pTrkTrkDtpos2,
                                               pTrkTrkDz,pTrkTrkDt);
              match = match && MatchOtherTracksInEvent(ac,vldcptr,track,objectlist);
            }
          }
          if(match){
            thiseventmatches = true;
            MSG("EventSR",Msg::kVerbose)
              << "  adding object to event \n";
            AddObjectToEvent(reco,objectlist,prevlist,merge);
            if(!merge) prevlist=objectlist;
            merge=true;
          }
        } // end loop over objects in event
      } // end loop over event list
        // we found no matches for this object, so start a new event
      if (!merge) {
        Bool_t pass=true;
        if (shower) {
          MSG("EventSR",Msg::kVerbose)
            << "    none found, creating event for shower\n";
        }
        
        if (track) {
          MSG("EventSR",Msg::kVerbose)
            << "    none found, creating event for track\n";
        }
        
        if(pass){
          TObjArray *event = new TObjArray;
          event->Add(recolist.At(i));
          eventlist.Add(event);
        }
      }
    } // end loop over recolist
  
    // 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);
      }
    }
    eventlist.Compress();
    
    // merge showers in NC events
    MergeShowers(cxx, eventlist,ac,slice,clusterlist,showerlist);
    
    // 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("EventSR",Msg::kVerbose)  << "making event of " 
                                    << objectlist->GetLast()+1 
                                    << " objects \n";
      cxx.SetDataIn(objectlist);
      buildEvent=true;
      CandEventHandle eventhandle = CandEvent::MakeCandidate(ah,cxx);
      MSG("EventSR",Msg::kVerbose) << " # of strips:  " 
                                   << eventhandle.GetNDaughters() << "\n";
      eventhandle.SetCandSlice(slice);
      ch.AddDaughterLink(eventhandle);
      delete eventlist.At(i);
    }
    if (!buildEvent) BuildEventFromUnassoc(ac,ch,cx,slice);
  }
  
  delete trackItr;
  delete showerItr;
  delete clusterItr;
    
  MSG("EventSR",Msg::kVerbose) << "starting unassociated hits assignment \n";

  // find unassociated hits, and place them in tobjarray
  sliceItr.Reset(); 
  TObjArray unassociated;
  CandShowerHandleItr * showerItr2=0;
  CandTrackHandleItr * trackItr2=0;
  CandEventListHandle &eventlist = dynamic_cast<CandEventListHandle &>(ch);
  CandEventHandleItr * eventItr2=0;
  if (showerlist && showerlist->GetNDaughters()>0) {
    showerItr2 = new CandShowerHandleItr(showerlist->GetDaughterIterator());
  }
  if (tracklist && tracklist->GetNDaughters()>0) { 
    trackItr2 = new CandTrackHandleItr(tracklist->GetDaughterIterator());
  }
  if (eventlist.GetNDaughters()>0) {    
    eventItr2 = new CandEventHandleItr(eventlist.GetDaughterIterator());
  }
  FindUnassociated(ac, sliceItr,eventItr2, unassociated);  
  delete trackItr2;
  delete showerItr2;
  delete eventItr2;

  // reset primary shower assignment for each event (is this needed?)
  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);
  
  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("EventSR",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("EventSR",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("EventSR",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
        if(closeststrip->GetPlane()>minShwPlane-pMaxNewShwLen &&
           closeststrip->GetPlane()<maxShwPlane+pMaxNewShwLen ){
          
          MSG("EventSR",Msg::kVerbose) << "adding strip to shower \n"; 
          AddStripToEvent(closestevent,shower,clusterlist,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,
                            clusterlist,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("EventSR",Msg::kVerbose) << "added " << n_found << " of "
                               << unassociated.GetLast()+1 
                               << " unassociated hits " << endl;

  cxx.SetDataIn(showerlist);  // Pass along persistable CandShowerList
  ReConstructShowers(ac,ch,cxx);

// For each event, we set it's primary shower to that which has the
// highest energy.

  SetPrimaryShowers(ch,ac);
}

//____________________________________________________________________
void AlgEventSRList::AddStripToEvent(CandEventHandle * closestevent,
                                     CandShowerHandle * shower,
                                     CandClusterListHandle * clusterlist, 
                                     CandStripHandle * closeststrip)
{
  closestevent->AddDaughterLink(*closeststrip);
  shower->AddDaughterLink(*closeststrip);
  Int_t foundcluster=0;
  for (Int_t icluster=0; icluster < shower->GetLastCluster()+1 &&
         foundcluster==0; icluster++) {
    CandClusterHandle *cl = const_cast<CandClusterHandle *>
      (shower->GetCluster(icluster));
    // add strip to cluster of the same view
    if (cl->GetPlaneView()==closeststrip->GetPlaneView()) {
      foundcluster=1;
      MSG("EventSR",Msg::kVerbose) << "adding strip to cluster "
                                   << closeststrip->GetPlane() << endl;
      TIter cchItr(clusterlist->GetDaughterIterator());
      while (CandClusterHandle *cch = 
             dynamic_cast<CandClusterHandle*>(cchItr())) {
        if(cch->IsEquivalent(cl)) {
          shower->RemoveCluster(cl);
          CandClusterHandle * newcl = cch->DupHandle();
          newcl->SetCandSlice(cch->GetCandSlice());
          newcl->AddDaughterLink(*closeststrip);
          clusterlist->RemoveDaughter(cch);
          clusterlist->AddDaughterLink(*newcl);
          shower->AddCluster(newcl);
          delete newcl;
          break;
        }
      }
    }
  }
}
//__________________________________________________________

void AlgEventSRList::CreatePrimaryShower(AlgConfig & ac, CandContext & cxx, 
                                         CandEventHandle * closestevent,
                                         CandShowerListHandle * showerlist,
                                         CandClusterListHandle * clusterlist,
                                         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
  
  // 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("EventSR",Msg::kError) <<"CandShowerListHandle missing\n" << endl;
    return;
  }
  
  closestevent->AddDaughterLink(*closeststrip);
  
  // create arrays of strips for each planeview to create
  // shower from
  
  MSG("EventSR",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
  if(ustriparray.GetLast()>=0 && vstriparray.GetLast()>=0){
    // now construct CandClusters
    AlgHandle clusterah = af.GetAlgHandle("AlgClusterSR","default");
    
    cxx.SetDataIn(&ustriparray);
    CandClusterHandle uclusterhandle =
      CandCluster::MakeCandidate(clusterah,cxx);
    uclusterhandle.SetCandSlice(closestevent->GetCandSlice());
    uclusterhandle.IsShowerLike(1);
    clusterlist->AddDaughterLink(uclusterhandle);
    
    cxx.SetDataIn(&vstriparray);
    CandClusterHandle vclusterhandle =
      CandCluster::MakeCandidate(clusterah,cxx);
    vclusterhandle.SetCandSlice(closestevent->GetCandSlice());
    vclusterhandle.IsShowerLike(1);
    clusterlist->AddDaughterLink(vclusterhandle);
    
    // build shower from U/V clusters just created
    TObjArray newshower;
    CandHandle * ucluster =
      clusterlist->FindDaughter(&uclusterhandle);
    CandHandle * vcluster =
      clusterlist->FindDaughter(&vclusterhandle);
    
    newshower.Add(ucluster);
    newshower.Add(vcluster);
    
    cxx.SetDataIn(&newshower);
    AlgHandle showerah = af.GetAlgHandle("AlgShowerSR",pEventAlgConfig);
    CandShowerHandle showerhandle =
      CandShower::MakeCandidate(showerah,cxx);
    showerhandle.SetCandSlice(closestevent->GetCandSlice());
    showerlist->AddDaughterLink(showerhandle);
    CandHandle *shw = showerlist->FindDaughter(&showerhandle);
    CandShowerHandle *shwh = dynamic_cast<CandShowerHandle*>(shw);
    
    // add this new shower to this event, and make it the
    // primary shower
    closestevent->AddShower(shwh);    
    closestevent->SetPrimaryShower(shwh);
  }
}

//____________________________________________________________________
void AlgEventSRList:: 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 AlgEventSRList::ReConstructShowers(AlgConfig &ac, CandHandle &ch, CandContext &cxx)
{

  // Retrieve persistant modifiable CandShowerList
  CandShowerListHandle *showerlist = const_cast<CandShowerListHandle*>
    (dynamic_cast<const CandShowerListHandle*>(cxx.GetDataIn()));

  // 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 ah = af.GetAlgHandle("AlgShowerSR",pEventAlgConfig);
  AlgConfig acshw = ah.GetAlgConfig();

  TIter evItr(ch.GetDaughterIterator());
  while (CandEventHandle *ev =
         dynamic_cast<CandEventHandle*>(evItr())) {
    for (Int_t ishower=0; ishower<ev->GetLastShower()+1; ishower++) {
      CandShowerHandle *shower =
        dynamic_cast<CandShowerHandle*>(ev->GetShowerWritable(ishower));
      if (shower) {
        TObjArray newshower;
        // fill  array of clusters for this shower
        for (Int_t icluster=0; icluster<shower->GetLastCluster()+1;
             icluster++) {
          const CandClusterHandle *cl = shower->GetCluster(icluster);
          if (cl) {
            // New handle on heap
            CandClusterHandle *cclh = cl->DupHandle();
            // Add new heap-based CandHandle to TObjArray
            newshower.Add(cclh);
          }
        }
        // re-run algorithm for this shower
        if (newshower.GetEntriesFast() > 0) {
          cxx.SetDataIn(&newshower);
          ah.RunAlg(*shower,cxx);
        }
        else {
          MSG("EventSR", Msg::kError)
            << "Attempt to pass empty TObjArray newshower to AlgShowerSR"
            << " in CandContext.  AlgShowerSR::RunAlg() call ignored."
            << endl;
        }

        // Delete heap-based CandHandles in TObjArray
        newshower.Delete();
      }
      else {
        MSG("EventSR",Msg::kWarning)
          << "Handle is not a CandShowerHandle."
          << "  Shower not processed in ReConstructShowers."
          << endl;
        continue;
      }

      // 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.
      
      CandTrackHandle *primaryTrack = ev->GetPrimaryTrack();
      
      // if event has primary track, extrapolate through shower and subtract energy loss to recalibrate 
      if (primaryTrack) shower->CalibrateEnergy(primaryTrack,acshw);
      
      // 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("EventSR",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);
        }
      }
    }
  }
}

//______________________________________________________________________
void AlgEventSRList::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 AlgEventSRList::BuildEventFromUnassoc(AlgConfig &ac,
                                           CandHandle &ch,
                                           CandContext &cx, 
                                           CandSliceHandle * slice)
{
  assert(cx.GetDataIn());
  MSG("EventSR",Msg::kDebug) <<" Building Event - no showers or tracks " 
                             << endl;
  if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
    MSG("EventSR",Msg::kDebug)<<" dataIn inherits from TObjArray " << endl;
    return;
  }
  
  const CandSliceListHandle *slicelist = 0;
  CandShowerListHandle *showerlist = 0;
  CandClusterListHandle *clusterlist = 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("CandClusterListHandle")) {
      clusterlist = dynamic_cast<CandClusterListHandle*>(tobj);
    }
  }
  
  if (!showerlist || !slicelist) {
    MSG("EventSR",Msg::kDebug)<<" no showerlist " << 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("EventSR",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("EventSR",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("EventSR",Msg::kDebug) << " 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("EventSR",Msg::kDebug)<< " 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("EventSR",Msg::kDebug) << " 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("EventSR",Msg::kDebug)<< " 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("EventSR",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("EventSR",Msg::kDebug) << " 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("EventSR",Msg::kDebug)<< " 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("EventSR",Msg::kDebug)<< " 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("EventSR",Msg::kDebug)<< " 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("EventSR",Msg::kDebug)<< " considering event formation  u/v: " 
                              << !firstu  << "/" << !firstv << " charge " 
                              << totcharge << endl;
    if (!firstu && !firstv && totcharge>20.0) {      
      found=true;
      
      Bool_t newShowerOverlapsOld=false;
      CandHandle *shw = NULL;

      // construct CandClusters
      cxx.SetDataIn(&neweventu);
      AlgHandle clusterah = af.GetAlgHandle("AlgClusterSR","default");
      CandClusterHandle uclusterhandle =
        CandCluster::MakeCandidate(clusterah,cxx);
      uclusterhandle.SetCandSlice(slice);
      uclusterhandle.IsShowerLike(1);
      clusterlist->AddDaughterLink(uclusterhandle);
      cxx.SetDataIn(&neweventv);
      CandClusterHandle vclusterhandle =
        CandCluster::MakeCandidate(clusterah,cxx);
      vclusterhandle.SetCandSlice(slice);
      vclusterhandle.IsShowerLike(1);
      clusterlist->AddDaughterLink(vclusterhandle);
      CandHandle *ucluster = clusterlist->FindDaughter(&uclusterhandle);
      CandHandle *vcluster = clusterlist->FindDaughter(&vclusterhandle);
      
      //build shower
      TObjArray newshower;
      newshower.Add(ucluster);
      newshower.Add(vcluster);
      cxx.SetDataIn(&newshower);
      
      
      AlgHandle showerah = af.GetAlgHandle("AlgShowerSR",pEventAlgConfig);
      CandShowerHandle showerhandle =
        CandShower::MakeCandidate(showerah,cxx);
      showerhandle.SetCandSlice(slice);
      showerhandle.SetCandRecord(slicelist->GetCandRecord());
      
      CandShowerHandleItr showerItr(showerlist->GetDaughterIterator());
      while (CandShowerHandle *shower =
             dynamic_cast<CandShowerHandle*>(showerItr())){
        if(showerhandle.BelongsWithShower(shower,ac,
                                          vldcptr,
                                          pShwShwDtpos2,
                                          pShwShwDz,
                                          pShwShwDt)){
          newShowerOverlapsOld=true;
          MSG("EventSR",Msg::kDebug) 
            << " new shower overlaps with previous - don't make new event " 
            << endl;
          break;
        }
      }
      
      if(!newShowerOverlapsOld){
        MSG("EventSR",Msg::kDebug)<< " Building Event! " << endl;
        showerlist->AddDaughterLink(showerhandle);
        shw = showerlist->FindDaughter(&showerhandle);
      }
      
      if(!newShowerOverlapsOld){
        // form event from shower
        TObjArray recolist;
        recolist.Add(shw);
        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);
      }
    }
  }
}


//_________________________________________________________________________
void AlgEventSRList::RemoveTracksinShowers(AlgConfig &ac, TObjArray * recolist){


  Double_t pMinTrackIsolation = ac.GetDouble("MinTrackIsolation");
  Double_t pMinTrackIsolationDist = ac.GetDouble("MinTrackIsolationDist");

  for (Int_t i = 0; i <= recolist->GetLast(); i++) {
    CandTrackHandle *track=0;
    Int_t nisolated = 0;
    Int_t nisolatedU = 0;
    Int_t nisolatedV = 0;
    Bool_t isolated = true;
    Int_t longestcontig = 0;
    Int_t longestcontigU = 0;
    Int_t longestcontigV = 0;
 
    if (recolist->At(i)->InheritsFrom("CandTrackHandle")) {
      track = dynamic_cast<CandTrackHandle*>(recolist->At(i));
      
      // count number of planes in each view in which
      // track is outside shower
      Int_t firstplane = TMath::Min(track->GetBegPlane(),
                                    track->GetEndPlane());
      Int_t lastplane = TMath::Max(track->GetBegPlane(),
                                   track->GetEndPlane());
      Int_t firstuplane = TMath::Min(track->GetBegPlane(PlaneView::kU),
                                     track->GetEndPlane(PlaneView::kU));
      Int_t firstvplane = TMath::Min(track->GetBegPlane(PlaneView::kV),
                                     track->GetEndPlane(PlaneView::kV));
      Int_t lastuplane = TMath::Max(track->GetBegPlane(PlaneView::kU),
                                  track->GetEndPlane(PlaneView::kU) );
      Int_t lastvplane = TMath::Max(track->GetBegPlane(PlaneView::kV),
                                  track->GetEndPlane(PlaneView::kV) );      

      MSG("RemoveTracksInShowers",Msg::kDebug) << " track extent " << firstplane << " " << lastplane <<  " " << track->GetNDaughters() << " " << track->GetBegPlane() << " " << track->GetEndPlane() << endl;

      for(Int_t iplane = firstplane;iplane <= lastplane;iplane++){
        PlexPlaneId plnid(track->GetVldContext()->GetDetector(),iplane,false);
        // require that the plane be within the track 
        // extent for this planeview
        isolated = true;
        if((plnid.GetPlaneView() == PlaneView::kU  && 
            iplane >= firstuplane && iplane <= lastuplane) ||
           (plnid.GetPlaneView() == PlaneView::kV  && 
            iplane >= firstvplane && 
            iplane<=lastvplane)){
        
          for (Int_t j = 0; j<=recolist->GetLast(); j++) {
            if(isolated){
              if (recolist->At(j)->InheritsFrom("CandShowerHandle")) {
                CandShowerHandle * shower = dynamic_cast<CandShowerHandle*>(recolist->At(j));
                if(iplane >= shower->GetBegPlane() &&
                   iplane <= shower->GetEndPlane()){
                  // if plane within shower extent, compare track position
                  // with shower width. If track outside shower by
                  // a distance pMinTrackIsolationDist, or shower width
                  // is less than this distance, inc isolated.
                  Double_t minPE=2.0;
                  Float_t minUshw = shower->GetMinU(iplane,minPE)+0.02;
                  Float_t minVshw = shower->GetMinV(iplane,minPE)+0.02;
                  Float_t maxUshw = shower->GetMaxU(iplane,minPE)-0.02;
                  Float_t maxVshw = shower->GetMaxV(iplane,minPE)-0.02;

                  MSG("RemoveTracksInShowers",Msg::kDebug) << " shower extent Plane: " << iplane << " U " << minUshw << "/" <<  track->GetU(iplane) << "/" << maxUshw << " V  " << minVshw << "/" << track->GetV(iplane) << "/" << maxVshw << endl;
                  
                  if(plnid.GetPlaneView() == PlaneView::kU && 
                     (track->GetU(iplane) >= (minUshw-pMinTrackIsolationDist) && 
                      track->GetU(iplane) <= (maxUshw+pMinTrackIsolationDist) &&
                      maxUshw-minUshw > (pMinTrackIsolationDist))){
                    isolated = false;
                  }
                  if(plnid.GetPlaneView() == PlaneView::kV &&
                     (track->GetV(iplane) >= (minVshw-pMinTrackIsolationDist)&& 
                      track->GetV(iplane) <= (maxVshw+pMinTrackIsolationDist)&&
                      maxVshw-minVshw > (pMinTrackIsolationDist))){
                    isolated = false;
                  }
                  if(!track->IsTPosValid(iplane))isolated = false;
                }
              }      
            }
          }
          if(track->IsTPosValid(iplane)){
            nisolated++;
            if(plnid.GetPlaneView() == PlaneView::kU){
              nisolatedU++;
            }
            else{
              nisolatedV++;
            }
          }
          if(!isolated) {
            nisolated=0;  // reset contig. isolation counter
            if(plnid.GetPlaneView() == PlaneView::kU){
              nisolatedU=0;
            }
            else{
              nisolatedV=0;
            }
          }
          if(nisolated > longestcontig)longestcontig = nisolated;
          if(nisolatedU > longestcontigU)longestcontigU = nisolatedU;
          if(nisolatedV > longestcontigV)longestcontigV = nisolatedV;
           
          MSG("RemoveTracksInShowers",Msg::kDebug) << " nisolated, longest " << nisolated << " " << longestcontig << endl;

          // if the number of contiguous isolated track planes is less than threshold
          // we remove this track from the event         
        }
      } 
    }
    MSG("RemoveTracksInShowers",Msg::kDebug) << " longest contig:" << longestcontig  << " U " << longestcontigU << " V " << longestcontigV<< "/" << pMinTrackIsolation << endl;
    if(track && longestcontig < pMinTrackIsolation &&
       longestcontigU < pMinTrackIsolation &&
       longestcontigV < pMinTrackIsolation){
      MSG("RemoveTracksInShowers",Msg::kDebug) << " REMOVING TRACK" << endl;
      recolist->RemoveAt(i);     
      recolist->Compress();
      //If you compress here then decrement i so that you 
      //don't miss any tracks in recolist
      i-=1;
    }
  }
}
//_________________________________________________________________________________
// for a given slice, fill an array with all tracks and showers
void AlgEventSRList::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("EventSR", 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("EventSR", Msg::kDebug)<< "Finder Track being used " << endl;
          }
          recolist.Add(track);
        }
      }
      else {
        MSG("EventSR", Msg::kError)
          << "Track doesn't contain a slice.  Not added to recolist."
          << endl;
      }
    }
  }
}

//_________________________________________________________________________
// compare pre-existing tracks in object list with  new track, 
// and return false if they are inconsistent with being in the
// same event.  Since these tracks are tentatively determined to be
// from the same event, loosen requirements by factor of 3. 
Bool_t AlgEventSRList::MatchOtherTracksInEvent(AlgConfig &ac, const VldContext * vldcptr, 
                                               CandTrackHandle * track, TObjArray * objectlist){

  Double_t pTrkTrkDtpos2 = ac.GetDouble("TrkTrkDtpos2")*9;
  Double_t pTrkTrkDz = ac.GetDouble("TrkTrkDz")*3;
  Double_t pTrkTrkDt = ac.GetDouble("TrkTrkDt");
  for (Int_t ireco2=0; ireco2<=objectlist->GetLast();
       ireco2++) {
    if (objectlist->At(ireco2)->InheritsFrom("CandTrackHandle")) {
      CandTrackHandle *othertrack = dynamic_cast<CandTrackHandle*>
        (objectlist->At(ireco2));
      if( !track->BelongsWithTrack(othertrack, 
                                   ac, 
                                   vldcptr, 
                                   pTrkTrkDtpos2,pTrkTrkDz,pTrkTrkDt))return false;
    } 
  }
  return true;
}

//_______________________________________________________________
void  AlgEventSRList::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("EventSR",Msg::kVerbose) << "    adding to list\n";
  }
  else {                  
    MSG("EventSR",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();
  }
}

//____________________________________________________________
// for events with only showers, merge all showers

void AlgEventSRList::MergeShowers( CandContext &cxx, TObjArray & eventlist, 
                                   AlgConfig &ac,
                                   CandSliceHandle * slice,
                                   CandClusterListHandle * clusterlist, 
                                   CandShowerListHandle * showerlist){
  
  
  const char *pEventAlgConfig = 0;
  ac.Get("EventAlgConfig",pEventAlgConfig);

  for (Int_t i=0; i<=eventlist.GetLast(); i++) {
    TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(i));
    if(objectlist->GetLast()<=0)break;
    Bool_t foundtrack=false;
    for (Int_t i=0; i<=objectlist->GetLast(); i++) {
      if (!(objectlist->At(i)->InheritsFrom("CandShowerHandle"))) {
        foundtrack=true;
        break;
      }
    }
    if(foundtrack)break;
    MSG("EventSR",Msg::kDebug) << "  merging showers\n";
    
    // if we are here, we have an event with multiple showers and no tracks, so we merge into single shower.
    
    // Get singleton instance of AlgFactory
    AlgFactory &af = AlgFactory::GetInstance();
    
    // construct total U/V cluster daughter list
    TObjArray ustriparray;
    TObjArray vstriparray;
    for (Int_t i=0; i<=objectlist->GetLast(); i++) {
      CandShowerHandle * shower = dynamic_cast<CandShowerHandle *>(objectlist->At(i));
      TIter stripItr(shower->GetDaughterIterator());
      while (CandStripHandle *strip =
             dynamic_cast<CandStripHandle*>(stripItr())) {
        if(strip->GetPlaneView()==PlaneView::kU)ustriparray.Add(strip);
        if(strip->GetPlaneView()==PlaneView::kV)vstriparray.Add(strip);
      }
    }
    
    // make U cluster
    cxx.SetDataIn(&ustriparray);
    AlgHandle clusterah = af.GetAlgHandle("AlgClusterSR","default");
    CandClusterHandle uclusterhandle =
      CandCluster::MakeCandidate(clusterah,cxx);
    uclusterhandle.SetCandSlice(slice);
    uclusterhandle.IsShowerLike(1);
    
    // make V cluster
    cxx.SetDataIn(&vstriparray);
    CandClusterHandle vclusterhandle =
      CandCluster::MakeCandidate(clusterah,cxx);
    vclusterhandle.SetCandSlice(slice);
    vclusterhandle.IsShowerLike(1);
    
    // add new clusters to cluster list
    clusterlist->AddDaughterLink(uclusterhandle);
    clusterlist->AddDaughterLink(vclusterhandle);
    
    // build shower
    TObjArray newshower;
    CandHandle *ucluster =
      clusterlist->FindDaughter(&uclusterhandle);
    CandHandle *vcluster =
      clusterlist->FindDaughter(&vclusterhandle);
    newshower.Add(ucluster);
    newshower.Add(vcluster);
    cxx.SetDataIn(&newshower);

    AlgHandle showerah = af.GetAlgHandle("AlgShowerSR",pEventAlgConfig);
    CandShowerHandle showerhandle =
      CandShower::MakeCandidate(showerah,cxx);
    showerhandle.SetCandSlice(slice);
    showerlist->AddDaughterLink(showerhandle);
    // need to clean up shower list
    for (Int_t i=0; i<=objectlist->GetLast(); i++) {
      CandShowerHandle * shower = dynamic_cast<CandShowerHandle *>(objectlist->At(i));
      TIter shiter(showerlist->GetDaughterIterator());
      CandShowerHandle *oldshower;
      Bool_t found=false;
      while (!found &&
             (oldshower = dynamic_cast<CandShowerHandle*>(shiter()))) {
        if (shower->IsEqual(oldshower)) {
          found=true;
          showerlist->RemoveDaughter(oldshower);
        }
      }  
    }    
    CandShowerHandleItr showerItr(showerlist->GetDaughterIterator());
    while (CandShowerHandle *shower =
           dynamic_cast<CandShowerHandle*>(showerItr())){
      MSG("EventSR",Msg::kDebug) << " new shower " << shower << " " << shower->GetVldContext() << endl;
    }
    // reset object list to contain new shower
    objectlist->Clear();
    CandShowerHandle * shower = dynamic_cast<CandShowerHandle *>(showerlist->FindDaughter(&showerhandle));
    assert (shower);
    objectlist->Add(shower);
  }
}


//__________________________________________________________________________
// fill TOjArray with strips not in any object daughter list in this snarl
void AlgEventSRList::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())) {

      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);
      }
    }
  }
}

//__________________________________________________________________________
// fill TOjArray with strips not in any object daughter list in this snarl
void AlgEventSRList::FindUnassociated(AlgConfig& ac, CandSliceHandleItr & sliceItr, CandShowerHandleItr * showerItr, CandTrackHandleItr* trackItr, 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())) {

      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 (showerItr) {
          showerItr->Reset();
          while (CandShowerHandle *shower =
                 dynamic_cast<CandShowerHandle*>((*showerItr)())) {
            if (shower->GetCandSlice()) {
              if (slice->IsCloneOf(*(shower->GetCandSlice()))) {
                if (shower->FindDaughter(strip)) {
                  found = true;
                  break;
                }
              }
            }
          }
        }
        if (trackItr && !found) {
          trackItr->Reset();
          while (CandTrackHandle *track =
                 dynamic_cast<CandTrackHandle*>((*trackItr)())) {
            if (track->GetCandSlice()) {
              if (slice->IsCloneOf(*(track->GetCandSlice()))) {
                if (track->FindDaughter(strip)) {
                  found = true;
                  break;
                }
              }
            }
          }
        }
      }
      // this strip is not in any event, add to unassociated list
      if (!found) {
        unassociated.Add(strip);
      }
    }
  }
}

//____________________________________________________________________

void AlgEventSRList::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("EventSR",Msg::kVerbose) << i << " " 
        << strip->GetPlane() << " " << strip->GetStrip() << " 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("EventSR",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("EventSR",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 AlgEventSRList::Trace(const char * /* c */) const
{
}

---