TridModelMaker.cxx

Go to the documentation of this file.

#include "TridModelMaker.h"

#include "TridModelStrip.h"
#include "TridModelIntersect.h"
#include "TridModelTrack.h"
#include "TridModelShower.h"
#include "TridModelRecoStrip.h"
#include "TridModelCrate.h"
#include "TridModelPmt.h"
#include "TridModelStrip.h"
#include "TridModelSlice.h"
#include "TridControl.h"

#include <TObject.h>
#include <TObjArray.h>
#include <TClonesArray.h>

#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "CandDigit/CandDigitListHandle.h"

#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRSlice.h"

#include "MessageService/MsgService.h"
#include "DataUtil/GetCandidate.h"
#include "DataUtil/GetVldContext.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Plex/PlexHandle.h"

#include "MessageService/MsgService.h"

#include <cmath>

CVSID("$Id: TridModelMaker.cxx,v 1.8 2006/12/01 20:25:49 rhatcher Exp $");

void TridModelMaker::Prepare(const MomNavigator* mom)
{
  // Explore Mom. See what there is to see.

  // Is there a CandRecord?
  fCandRecord = dynamic_cast<const CandRecord*> 
    (mom->GetFragment("CandRecord","PrimaryCandidateRecord"));
  
  // Is there an Ntuple record?
  fNtpStRecord = dynamic_cast<const NtpStRecord*> 
    (mom->GetFragment("NtpStRecord","Primary"));
   
  // Find VldContext.
  std::vector<VldContext> contexts = DataUtil::GetVldContext(mom);
  if(contexts.size()==0) {
    MSG("TriD",Msg::kError) << "No data in MOM: Unable to find VldContext." << std::endl;
  } else {
    fContext = DataUtil::GetVldContext(mom)[0];
  }
}

Int_t TridModelMaker::CreateTrackModels(const MomNavigator* mom,
                                        TridModelList&      models)
{
  // Create track models.
  Int_t nmodel=0;

  const CandTrackListHandle* tracklist =
    DataUtil::GetCandidate<CandTrackListHandle>(mom,
                                                0,"CandTrackSRList");
  if (tracklist) {
    TIter trackItr(tracklist->GetDaughterIterator());
    while(CandTrackHandle* track = dynamic_cast<CandTrackHandle*>(trackItr())) {
      //cout << "Got a track! " << endl;
      if(track) {
        TridModelTrack* model = new TridModelTrack();
        model->AddCandidate(*track);
        for(int p=0;p<500;p++) {
          if(track->IsTPosValid(p)) {
            if(p<model->fFirstPlane) model->fFirstPlane = p;
            if(p>model->fLastPlane) model->fLastPlane = p;
            model->fN[p] = 1;
            model->fU[p] = track->GetU(p);
            model->fV[p] = track->GetV(p);
            model->fT[p] = track->GetT(p);          
          }
        }
        
        // Create description:
        model->fDescription = "Track: \n";
        model->fDescription+= Form("Distance: %f.1 m \n",track->GetdS());
        model->fDescription+= Form("Range:    %f.1 g/cm^2\n",track->GetRange());
        model->fDescription+= Form("Momentum: %f.1 GeV\n",track->GetMomentum());
        model->fDescription+= Form("Vertex: uvz=(%f.1, %f.1, %f.1)\n",
                                   track->GetVtxU(),track->GetVtxV(),track->GetVtxZ());
        model->fDescription+= Form("VtxDir: uvz=(%f.2, %f.2, %f.2)\n",
                                   track->GetVtxDirCosU(),track->GetVtxDirCosV(),track->GetVtxDirCosZ());
        model->fDescription+= Form("Dir:    uvz=(%f.2, %f.2, %f.2)\n",
                                   track->GetDirCosU(),track->GetDirCosV(),track->GetDirCosZ());
        model->fDescription+= Form("Planes: %d to %d\n",model->fFirstPlane,model->fLastPlane);
        model->fDescription+= Form("VtxDir: uvz=(%f.2,%f.2,%f.2)\n",
                                   track->GetVtxDirCosU(),track->GetVtxDirCosV(),track->GetVtxDirCosZ());
        model->fDescription+= Form("Vertex Trace: %.1f   TraceZ:%.1f \n",
                                   track->GetVtxTrace(), track->GetVtxTraceZ());
        model->fDescription+= Form("End Trace:    %.1f   TraceZ:%.1f \n",
                                   track->GetEndTrace(), track->GetEndTraceZ());
 
        models.AddModel(model);
        nmodel++;
      }
    }
  }

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->trk) {
      TIter iter(fNtpStRecord->trk);
       TObject* tobj;
       const NtpSRTrack* track;
       while( (tobj = iter.Next()) ) {
         if( (track = dynamic_cast<const NtpSRTrack*>(tobj)) ) {
           TridModelTrack* tmodel = new TridModelTrack;    
           for(int i=0;i<track->nstrip;i++) {
             int stripid = track->stp[i];
             const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(fNtpStRecord->stp->At(stripid));
             tmodel->AddStrip(stripid,
                              fContext,
                              strip,
                              (fContext.GetDetector()==Detector::kNear)?StripEnd::kWest:StripEnd::kWhole
                             );
             if(strip) {
               int plane = strip->plane;
               if(plane < tmodel->fFirstPlane) tmodel->fFirstPlane = plane;
               if(plane > tmodel->fLastPlane)  tmodel->fLastPlane  = plane;
               tmodel->fN[plane] += 1;
               tmodel->fU[plane] += track->stpu[i];
               tmodel->fV[plane] += track->stpv[i];
               double t = 0;
               double ends = 0;
               if(strip->time0>-999990) { t+=(strip->time0 - fContext.GetTimeStamp().GetNanoSec()*1e-9); ends++; }
               if(strip->time1>-999990) { t+=(strip->time1 - fContext.GetTimeStamp().GetNanoSec()*1e-9); ends++; }
               tmodel->fT[strip->plane] += t/ends;
               tmodel->fTotalCharge[CalStripType::kSigMapped] += track->stpph0sigmap[i];
               tmodel->fTotalCharge[CalStripType::kSigMapped] += track->stpph1sigmap[i];
               tmodel->fTotalCharge[CalStripType::kMIP] += track->stpph0mip[i];
               tmodel->fTotalCharge[CalStripType::kMIP] += track->stpph1mip[i];
               tmodel->fTotalCharge[CalStripType::kGeV] += track->stpph0gev[i];
               tmodel->fTotalCharge[CalStripType::kGeV] += track->stpph1gev[i];
             }
           }
           tmodel->fTotalCharge[CalStripType::kNone]      = track->ph.raw;
           tmodel->fTotalCharge[CalStripType::kSigLin]    = track->ph.siglin;
           tmodel->fTotalCharge[CalStripType::kSigCorr]   = track->ph.sigcor;
           tmodel->fTotalCharge[CalStripType::kPE]        = track->ph.pe;
           tmodel->fTotalCharge[CalStripType::kSigMapped] = track->ph.sigmap;
           tmodel->fTotalCharge[CalStripType::kMIP]       = track->ph.mip;
           tmodel->fTotalCharge[CalStripType::kGeV]       = track->ph.gev;

           // Normalize to number of strips per plane.
           for(int i=0;i<500;i++) 
             if(tmodel->fN[i]>0) {
               tmodel->fU[i]/=(float)tmodel->fN[i];
               tmodel->fV[i]/=(float)tmodel->fN[i];            
               tmodel->fT[i]/=(double)tmodel->fN[i];           
             }
         
           
           // Create description:
           tmodel->fDescription = "Track: \n";
           tmodel->fDescription+= Form("Distance: %f.1 m \n",track->ds);
           tmodel->fDescription+= Form("Range:    %f.1 g/cm^2\n",track->range);
           tmodel->fDescription+= Form("Momentum: %f.1 GeV (range)\n",track->range);
           tmodel->fDescription+= Form("q/p:      %f.1 +/- %.1f \n",track->momentum.qp, track->momentum.eqp);
           tmodel->fDescription+= Form("Vertex: uvz=(%f.1, %f.1, %f.1)\n",
                                       track->vtx.u, track->vtx.v, track->vtx.z );
           tmodel->fDescription+= Form("VtxDir: uvz=(%f.2, %f.2, %f.2)\n",
                                       track->vtx.dcosu, track->vtx.dcosv, track->vtx.dcosz);
           tmodel->fDescription+= Form("Dir:    uvz=(%f.2, %f.2, %f.2)\n",
                                       track->lin.dcosu, track->lin.dcosv, track->lin.dcosz);
           tmodel->fDescription+= Form("Planes: %d to %d\n",
                                       tmodel->fFirstPlane,tmodel->fLastPlane);
           tmodel->fDescription+= Form("Vertex Trace: %.1f   TraceZ:%.1f \n",
                                       track->fidvtx.trace, track->fidvtx.tracez);
           tmodel->fDescription+= Form("End Trace:    %.1f   TraceZ:%.1f \n",
                                       track->fidend.trace, track->fidend.tracez);
 
           models.AddModel(tmodel);
           nmodel++;
         }
       }
    }
  }


  return nmodel;
}


Int_t TridModelMaker::CreateShowerModels(const MomNavigator* mom,
                                         TridModelList& models)  
{
  // Create shower models.
  Int_t nmodel=0;
  
  UgliGeomHandle ugli(fContext);
  
  const CandShowerListHandle* showerlist =
    DataUtil::GetCandidate<CandShowerListHandle>(mom,
                                                 0,"CandShowerList");
  if (showerlist) {
    TIter showerItr(showerlist->GetDaughterIterator());
    while(CandShowerHandle* shower = dynamic_cast<CandShowerHandle*>(showerItr())) {
      //cout << "Got a shower! " << endl;
      if(shower) {
        TridModelShower* showermodel = new TridModelShower;
        showermodel->AddCandidate(*shower);
        for(int i=0;i<500;i++) {
          if(shower->GetNStrips(i)>0) {
            PlexPlaneId plane(fContext.GetDetector(),i);
            showermodel->fN[i] = shower->GetNStrips(i);
            showermodel->fU[i] = shower->GetU(i);
            showermodel->fV[i] = shower->GetV(i);

            if(plane.GetPlaneView()==PlaneView::kU) {
              showermodel->fWidth[i] = TMath::Max(fabs(shower->GetMinU(i)-shower->GetU(i)),
                                                  fabs(shower->GetMaxU(i)-shower->GetU(i)));
            } else {
              showermodel->fWidth[i] = TMath::Max(fabs(shower->GetMinV(i)-shower->GetV(i)),
                                                  fabs(shower->GetMaxV(i)-shower->GetV(i)));        }
          }
          showermodel->fT[i] = shower->GetT(i);
        }      
        models.AddModel(showermodel);
        nmodel++;

        // Create shower recostrips
        TIter stripItr(shower->GetDaughterIterator());
        while(CandStripHandle* strip = dynamic_cast<CandStripHandle*>(stripItr())) {
          UgliStripHandle ustrip = ugli.GetStripHandle(strip->GetStripEndId());
          Int_t plane = strip->GetPlane();
          Float_t width_u = 0.1;
          Float_t width_v = 0.1;
          Float_t nu = 0;
          Float_t nv = 0;
          for(int p = plane-1; p<= plane+1; p++) {
            if(shower->IsTPosValid(p)) {
              PlexPlaneId planeId(fContext.GetDetector(),p);
              if(planeId.GetPlaneView()==PlaneView::kU) {
                width_u += (shower->GetMaxU(p)-shower->GetMinU(p))*0.5; 
                nu++;
              } else {
                width_v += (shower->GetMaxV(p)-shower->GetMinV(p))*0.5;
                nv++;
              }
            }
          }
          if(nu>0) width_u /= nu;
          if(nv>0) width_v /= nv;

          TridModelRecoStrip* recostrip = 
            new TridModelRecoStrip(ustrip,
                                   shower->GetU(plane),
                                   shower->GetV(plane),
                                   width_u,
                                   width_v,
                                   TridModelRecoStrip::kShower,
                                   showermodel
                                   );
          recostrip->AddCandidate(*strip);

          recostrip->fTotalCharge[CalStripType::kSigMapped] += shower->GetStripCharge(strip,CalStripType::kSigMapped);
          recostrip->fTotalCharge[CalStripType::kMIP]    += shower->GetStripCharge(strip,CalStripType::kMIP);
          recostrip->fTotalCharge[CalStripType::kGeV]    += shower->GetStripCharge(strip,CalStripType::kGeV);

          models.AddModel(recostrip);
          nmodel++;
        }
      }
    }
  }

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->shw) {
      TIter iter(fNtpStRecord->shw);
       TObject* tobj;
       const NtpSRShower* shower;
       while( (tobj = iter.Next()) ) {
         if( (shower = dynamic_cast<const NtpSRShower*>(tobj)) ) {
           //if(shower->stpu == 0) { cout << "Bad shower in ntp! " << endl; continue; };
           TridModelShower* model = new TridModelShower;           
           for(int i=0;i<shower->nstrip;i++) {
             int stripid = shower->stp[i];
             const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(fNtpStRecord->stp->At(stripid));
             model->AddStrip(stripid,
                             fContext,
                             strip,
                             (fContext.GetDetector()==Detector::kNear)?StripEnd::kWest:StripEnd::kWhole
                             );
             if(strip) {
               PlexPlaneId plane(fContext.GetDetector(),strip->plane);
               model->fN[strip->plane] += 1;
               model->fU[strip->plane] += strip->tpos; // just use U as tpos holder
               double t = 0;
               double ends = 0;
               if(strip->time0>-999990) { t+=(strip->time0 - fContext.GetTimeStamp().GetNanoSec()*1e-9); ends++; }
               if(strip->time1>-999990) { t=+(strip->time1 - fContext.GetTimeStamp().GetNanoSec()*1e-9); ends++; }
               model->fT[strip->plane] += t/ends;
             }
             model->fTotalCharge[CalStripType::kNone]      = shower->ph.raw;
             model->fTotalCharge[CalStripType::kSigLin]    = shower->ph.siglin;
             model->fTotalCharge[CalStripType::kSigCorr]   = shower->ph.sigcor;
             model->fTotalCharge[CalStripType::kPE]        = shower->ph.pe;
             model->fTotalCharge[CalStripType::kSigMapped] = shower->ph.sigmap;
             model->fTotalCharge[CalStripType::kMIP]       = shower->ph.mip;
             model->fTotalCharge[CalStripType::kGeV]       = shower->ph.gev;
           }
         
           // Normalize to number of strips per plane.
           for(int p=0;p<500;p++) {
             if(model->fN[p]>0) {
               model->fU[p]/=(float)model->fN[p];
               if(model->fFirstPlane>p) model->fFirstPlane=p;
               if(model->fLastPlane<p)  model->fLastPlane=p;
               model->fT[p]/=(double)(model->fN[p]);
             }
           }

           // Use fV to store lpos.
           // Use last/next plane average in the middle of shower
           // use next or last plane at ends.
           for(int p=model->fFirstPlane+1; p<model->fLastPlane;p++) {
             model->fV[p] = 0.5* model->fU[p-1] + 0.5 * model->fU[p+1];      
           }
           model->fV[model->fFirstPlane] = model->fU[model->fFirstPlane+1];
           model->fV[model->fLastPlane]  = model->fU[model->fLastPlane-1];

           // Get the widths.
           for(int i=0;i<shower->nstrip;i++) {
             int stripid = shower->stp[i];
             const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(fNtpStRecord->stp->At(stripid));
             if(strip){
               PlexPlaneId plane(fContext.GetDetector(),strip->plane);
               float w;
               float center = model->fU[strip->plane];
               w = fabs(strip->tpos - center);
               if(w > model->fWidth[strip->plane]) model->fWidth[strip->plane] = w;
             }  
           }
           
           // Rescramble to get UV correct:      
           for(int p=model->fFirstPlane; p<model->fLastPlane+1;p++) {  
             PlexPlaneId plane(fContext.GetDetector(),p);
             if(plane.GetPlaneView()==PlaneView::kV){
               float u = model->fV[p];
               model->fV[p] = model->fU[p];
               model->fU[p] = u;
             }
             //cout << Form("%3d: %s %4.1f %4.1f %4.1f",
             //                   p,
             //           PlaneView::AsString(plane.GetPlaneView()),
             //           model->GetU(p),
             //           model->GetV(p),
             //           model->GetWidth(p)) << endl;
           }


           models.AddModel(model);
           nmodel++;
           
           // Create the shower recostrips.
           for(int i=0;i<shower->nstrip;i++) {
             int stripid = shower->stp[i];
             const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(fNtpStRecord->stp->At(stripid));
             if(strip){
               PlexStripEndId seid(fContext.GetDetector(),strip->plane,strip->strip);
               UgliStripHandle ustrip = ugli.GetStripHandle(seid);
               Int_t plane = strip->plane;
               float width_u = 0.04;
               float width_v = 0.04;
               Float_t nu = 0;
               Float_t nv = 0;
               for(int p = plane-1; p<= plane+1; p++) {
                 if(model->IsValid(p)) {
                   PlexPlaneId planeId(fContext.GetDetector(),p);
                   if(planeId.GetPlaneView()==PlaneView::kU) {
                     width_u += model->GetWidth(p);
                     nu++;
                   } else {
                     width_v += model->GetWidth(p);
                     nv++;
                   }
                 }
               }
               if(nu>0) width_u /= nu;
               if(nv>0) width_v /= nv;

               TridModelRecoStrip* recostrip = 
                 new TridModelRecoStrip(ustrip,
                                        model->GetU(strip->plane),
                                        model->GetV(strip->plane),
                                        width_u,
                                        width_v,
                                        TridModelRecoStrip::kShower,
                                        model
                                        );
               
               recostrip->AddStrip(strip->index,fContext,
                                   strip,
                                   (fContext.GetDetector()==Detector::kNear)?StripEnd::kWest:StripEnd::kWhole
                                   );
               models.AddModel(recostrip);
               nmodel++;
               
             }  
           }
           
         }
       }
    }
  }
  
  return nmodel;
}


Int_t TridModelMaker::CreateStripModels(const MomNavigator* mom,
                                        TridModelList& models)
{
  Int_t nmodel=0;
  PlexHandle plex(fContext);
  bool fBestDemuxOnly = true; // FIXME: configurable.
  
  const CandDigitListHandle* digitList =
    DataUtil::GetCandidate<CandDigitListHandle>(mom,
                                                "CandDigitListHandle",
                                                "canddigitlist");
 
  // Create strip models.

  if(digitList) {
    CandDigitHandleItr cdhItr(digitList->GetDaughterIterator());
    while(CandDigitHandle *cdh = cdhItr()) {
      PlexSEIdAltL altList = cdh->GetPlexSEIdAltL();
      
      int numitems = altList.GetSize();
      
    // Do only one strip
    // unless we're explicitly doing all
    // or these digit(s) have been selected.
      if((numitems>1)&&(fBestDemuxOnly))
        numitems = 1;
      
      for(int i=0; i<numitems; i++) {
        PlexSEIdAltLItem item = altList[i];
        if(numitems==1) item = altList.GetBestItem();
        
        // Make a new model, which might be redundant.
        TridModelStrip* newmodel = new TridModelStrip(item.GetSEId());
        
        // Try to find a matching one.
        bool need_new = true;
        TridModel* oldmodel;
        TridModelItr itr = models.GetIterator(newmodel->GetSortKey());
        while( (oldmodel = itr.Next()) ) {
          TridModelStrip* oldstrip = dynamic_cast<TridModelStrip*>(oldmodel);
          if(oldstrip) {
            if(oldstrip->ShouldContain(item.GetSEId()) ) {
              oldstrip->AddDigit(*cdh);
              need_new = false;
              delete newmodel;
              break;
            }
          }
        }
        
        if(need_new){
          newmodel->AddDigit(*cdh);
          models.AddModel(newmodel);    
          nmodel++;
          MSG("TriD",Msg::kDebug) << "Created strip model - "
                                  << " ModelAdd:"     << newmodel 
                                  << " Model ID:"     << ((newmodel) ? (newmodel->GetId()) : 0)
                                  << " ModelKey:"     << ((newmodel) ? (newmodel->GetSortKey()) : 0)
                                  << std::endl;
        }
      }
      
    }
  } // have digitList

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->stp) {
       TIter iter(fNtpStRecord->stp);
       TObject* tobj;
       const NtpSRStrip* strip;
       while( (tobj = iter.Next()) ) {
         if( (strip = dynamic_cast<const NtpSRStrip*>(tobj)) ) {
           PlexStripEndId seid(fContext.GetDetector(),
                               strip->plane,
                               strip->strip,
                               StripEnd::kUnknown);        

           if(strip->ph0.raw>0) {
             seid.SetEnd(StripEnd::kNegative);
             if(strip->ph1.raw>0) 
               seid.SetEnd(StripEnd::kWhole);
           } else if (strip->ph1.raw>0) 
             seid.SetEnd(StripEnd::kPositive);

           TridModelStrip* model = new TridModelStrip(seid);
           model->AddStrip(strip->index,fContext,strip,
                             (fContext.GetDetector()==Detector::kNear)?StripEnd::kWest:StripEnd::kWhole
                             );
           models.AddModel(model);
           nmodel++;
         }
       }
    }
  }

  return nmodel;
}


Int_t TridModelMaker::CreateIntersectionModels(const MomNavigator* /*mom*/,
                                               TridModelList& models)
{
  // Create intesection models
  
  Int_t nmodel = 0;
  
  // loop through existing models. Models should be sorted by plane number.
  TridModel*      model1;
  TridModel*      model2;
  TridModelStrip* modelStrip1;
  TridModelStrip* modelStrip2;

  TridModelItr itr = models.GetIterator();
  while( (model1 = itr.Next()) ) {
    modelStrip1 = dynamic_cast<TridModelStrip*>(model1);
    if(modelStrip1) {
      // Get the plane numbers.
      Int_t plane1 = modelStrip1->GetSortKey();
      Int_t plane2 = plane1+1;
      if(plane1<500) { // It's not in the veto shield

        // Ask for any models matching plane 2
        TridModelItr matchItr = models.GetIterator(plane2);
        while( (model2 = matchItr.Next() ) ) {
          modelStrip2 = dynamic_cast<TridModelStrip*>(model2);
          if(modelStrip2) {
            // We now have two strips in neigboring planes.
            // Are they in time?
            if(fabs(modelStrip1->GetMeanTime() - modelStrip2->GetMeanTime())<100e-9) {
              // Make the intersection.
              TridModelIntersect* intersectModel =
                new TridModelIntersect( *modelStrip1, *modelStrip2 );
              models.AddModel(intersectModel);
              nmodel++;
              MSG("TriD",Msg::kDebug) << "Created intersection model - "
                                      << " ModelAdd:"     << intersectModel 
                                      << " Model ID:"     << ((intersectModel) ? (intersectModel->GetId()) : 0)
                                      << " ModelKey:"     << ((intersectModel) ? (intersectModel->GetSortKey()) : 0)
                                      << std::endl;
              
              // Mark strips as having been correctly intersected.
              modelStrip1->fIntersections++;
              modelStrip2->fIntersections++;
            }
          }
        }
      }
    }
  }

  return nmodel;
}


Int_t TridModelMaker::CreatePmtModels(const MomNavigator* mom,
                                      TridModelList& models)
{
  // Build PMT list.
  Int_t nmodel=0;

  const CandDigitListHandle* digitList =
    DataUtil::GetCandidate<CandDigitListHandle>(mom,
                                                "CandDigitListHandle",
                                                "canddigitlist");

  if(digitList) {
    CandDigitHandleItr cdhItr(digitList->GetDaughterIterator());
    while(CandDigitHandle *cdh = cdhItr()) {
      PlexSEIdAltL altList = cdh->GetPlexSEIdAltL();
      if(altList.GetSize()>0) {
        PlexSEIdAltLItem item = altList.GetBestItem();
        
        // Make a model, or find the correct one.
        bool need_new_pmt = true;
        bool need_new_pixel = true;
        
        // Make a model to get the sort key.
        TridModelPmt newmodel(item.GetPixelSpotId());
        
        TridModel* oldmodel;
        TridModelItr itr = models.GetIterator(newmodel.GetSortKey());
        while( (oldmodel = itr.Next()) ) {
          TridModelPmt* oldpmt = dynamic_cast<TridModelPmt*>(oldmodel);
          if(oldpmt) 
            if(oldpmt->ShouldContain(item.GetPixelSpotId()) ) {
              // We already have this one. Add in this item.
              oldpmt->AddDigit(*cdh);
              need_new_pmt = false;
            }
          TridModelPixel* oldpix = dynamic_cast<TridModelPixel*>(oldmodel);
          if(oldpix)    
            if(oldpix->ShouldContain(item.GetPixelSpotId()) ) {
              // We already have this one. Add in this item.
              oldpix->AddDigit(*cdh);
              need_new_pixel = false;
            }  
        }

        if(need_new_pmt) {      
          TridModel* m = new TridModelPmt(item.GetPixelSpotId());
          m->AddDigit(*cdh);
          models.AddModel(m);     
          nmodel++;
        }
        if(need_new_pixel) {    
          TridModel* m = new TridModelPixel(item.GetPixelSpotId());
          m->AddDigit(*cdh);
          models.AddModel(m);
          nmodel++;
        }
        
      }
      
    }
  } // have digitList

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->stp) {
      PlexHandle plex(fContext);
      
      TIter iter(fNtpStRecord->stp);
      TObject* tobj;
      const NtpSRStrip* strip;
      while( (tobj = iter.Next()) ) {
        if( (strip = dynamic_cast<const NtpSRStrip*>(tobj)) ) {
          PlexStripEndId seid(fContext.GetDetector(),
                              strip->plane,
                              strip->strip,
                              StripEnd::kUnknown);         

          for(int iend=1; iend<=2; iend++) {
            float raw = strip->ph0.raw;
            StripEnd::StripEnd_t end = StripEnd::kNegative;
            if(iend==2) {
              end = StripEnd::kPositive;
              raw = strip->ph1.raw;
            }
            if(raw>0) {
              // Got an end.
              seid.SetEnd(end);
              PlexPixelSpotId psid = plex.GetPixelSpotId(seid);
               
              // Make a model, or find the correct one.
              bool need_new_pmt = true;
              bool need_new_pixel = true;
               
              // Make a model to get the sort key.
              TridModelPmt newmodel(psid);
               
              TridModel* oldmodel;
              TridModelItr itr = models.GetIterator(newmodel.GetSortKey());
              while( (oldmodel = itr.Next()) ) {
                TridModelPmt* oldpmt = dynamic_cast<TridModelPmt*>(oldmodel);
                if(oldpmt) 
                  if(oldpmt->ShouldContain(psid) ) {
                    // We already have this one. Add in this item.
                    oldpmt->AddStrip(strip->index, fContext, strip, end);
                    need_new_pmt = false;
                  }
                TridModelPixel* oldpix = dynamic_cast<TridModelPixel*>(oldmodel);
                if(oldpix)      
                  if(oldpix->ShouldContain(psid) ) {
                    // We already have this one. Add in this item.
                    oldpix->AddStrip(strip->index, fContext, strip, end);
                    need_new_pixel = false;
                  }  
              }
        
              if(need_new_pmt) {        
                TridModelPmt* newpmt = new TridModelPmt(psid);
                newpmt->AddStrip(strip->index, fContext, strip, end);
                models.AddModel(newpmt);         
                nmodel++;
              }
              if(need_new_pixel) {      
                TridModelPixel* newpix = new TridModelPixel(psid);
                newpix->AddStrip(strip->index, fContext, strip, end);
                models.AddModel(newpix);
                nmodel++;
              }        
            }
          }
        }
      }
    }
  }
  return nmodel;
}

Int_t TridModelMaker::CreateChannelModels(const MomNavigator* mom,
                                          TridModelList& models)
{
  Int_t nmodel=0;

  const CandDigitListHandle* digitList =
    DataUtil::GetCandidate<CandDigitListHandle>(mom,
                                                "CandDigitListHandle",
                                                "canddigitlist");

  if(digitList) {

  CandDigitHandleItr cdhItr(digitList->GetDaughterIterator());
  while(CandDigitHandle *cdh = cdhItr()) {
    PlexSEIdAltL altList = cdh->GetPlexSEIdAltL();
    if(altList.GetSize()>0) {
      PlexSEIdAltLItem item = altList.GetBestItem();
      
      // Make a model, or find the correct one.
      
      // Make a new model, which might be redundant.
      TridModelCrate* newmodel = new TridModelCrate(cdh->GetChannelId());
      
      // Try to find a matching one.
      bool need_new = true;
      TridModelCrate* oldmodel;
      TridModelItr itr = models.GetIterator(newmodel->GetSortKey());
      while( (oldmodel = dynamic_cast<TridModelCrate*>(itr.Next())) ) {
        if(oldmodel->ShouldContain(cdh->GetChannelId()) ) {
          oldmodel->AddDigit(*cdh);
          need_new = false;
          delete newmodel;
          break;
        }
      }
        
      if(need_new) {
        newmodel->AddDigit(*cdh);
        models.AddModel(newmodel);      
        nmodel++;
      }
    }
  }
  }

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->stp) {
      PlexHandle plex(fContext);
      
      TIter iter(fNtpStRecord->stp);
      TObject* tobj;
      const NtpSRStrip* strip;
      while( (tobj = iter.Next()) ) {
        if( (strip = dynamic_cast<const NtpSRStrip*>(tobj)) ) {
          PlexStripEndId seid(fContext.GetDetector(),
                              strip->plane,
                              strip->strip,
                              StripEnd::kUnknown);         

          for(int iend=1; iend<=2; iend++) {
            float raw = strip->ph0.raw;
            StripEnd::StripEnd_t end = StripEnd::kNegative;
            if(iend==2) {
              end = StripEnd::kPositive;
              raw = strip->ph1.raw;
            }
            if(raw>0) {
              // Got an end.
              seid.SetEnd(end);
              RawChannelId rcid = plex.GetRawChannelId(seid);
              
              // Make a new model, which might be redundant.
              TridModelCrate* newmodel = new TridModelCrate(rcid);            
              
              // Try to find a matching one.
              bool need_new = true;
              TridModelCrate* oldmodel;
              TridModelItr itr = models.GetIterator(newmodel->GetSortKey());
              while( (oldmodel = dynamic_cast<TridModelCrate*>(itr.Next())) ) {
                if(oldmodel->ShouldContain(rcid) ) {
                  oldmodel->AddStrip(strip->index,fContext,
                                     strip, end);
                  need_new = false;
                  delete newmodel;
                  break;
                }
              }
              
              if(need_new) {
                newmodel->AddStrip(strip->index,fContext,
                                   strip, end);
                models.AddModel(newmodel);      
                nmodel++;
              }
            }
            
          }
        }
      }
    }
  }
  return nmodel;
}


Int_t TridModelMaker::CreateSliceModels(const MomNavigator* /* mom*/,
                                        TridModelList& models)
{
  Int_t nmodel=0;
  Int_t maxplane = 120;
  if(fContext.GetDetector()==Detector::kFar) maxplane = 300;

  // Try from ntuple.
  if((nmodel==0) && (fNtpStRecord)) {
    if(fNtpStRecord->slc) {
      TIter iter(fNtpStRecord->slc);
      TObject* tobj;
      const NtpSRSlice* slice;
      int islice = 0;
      while( (tobj = iter.Next()) ) {
        if( (slice = dynamic_cast<const NtpSRSlice*>(tobj)) ) {
          TridModelSlice* model = new TridModelSlice;
          model->fSlice = islice++;
          for(int i=0;i<slice->nstrip;i++) {
             int stripid = slice->stp[i];
             const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(fNtpStRecord->stp->At(stripid));
             model->AddStrip(stripid, fContext, strip,
                             (fContext.GetDetector()==Detector::kNear)?StripEnd::kWest:StripEnd::kWhole
                             );
             if(strip){        
               if(strip->plane < maxplane) {
                 model->fHits.push_back(
                                        TridModelSlice::SliceHit(strip->plane,strip->time1-fContext.GetTimeStamp().GetNanoSec()*1e-9));
               }
             }
          }             
          model->SetSortKey((int)(model->GetMeanTime()*10000000)); // To sort on time to the nearest 0.1 us
          models.AddModel(model);
          nmodel++;       
        }
      }
    }
  }
  return nmodel;
}

---