UberModuleLite.cxx

Go to the documentation of this file.

//
// 29.05.2003 -- MAK: Changed NtpSRFitVertex to NtpSRVertex
// uggg!  This doesn't work with R0.18.0
//

#include <iomanip>
#include <set>
#include "TROOT.h"
#include "TObjectTable.h"
#include "MINF_Classes/MINFast.h"
#include "CandData/CandRecord.h"
#include "CandData/CandHeader.h"
#include "CalDetSI/CandCalDetSIHandle.h"
#include "CandDigit/CandDigitListHandle.h"
#include "CandDigit/CandDigitHandle.h"
#include "REROOT_Classes/REROOT_Event.h"
#include "REROOT_Classes/REROOT_FLSDigit.h"
#include "REROOT_Classes/REROOT_StdHep.h"
#include "RerootExodus/RerootExodus.h"
#include "Plex/PlexStripEndId.h"
#include "Plex/PlexSEIdAltLItem.h"
#include "Plex/PlexHandle.h"
#include "Plex/PlexPixelSpotId.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Navigation/NavSet.h"
#include "Conventions/CalDigitType.h"
#include "Conventions/CalTimeType.h"
#include "Conventions/Munits.h"
#include "Conventions/ReadoutType.h"
#include "Calibrator/Calibrator.h"
#include "MinosObjectMap/MomNavigator.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "JobControl/JobCResult.h"
#include "JobControl/JobCModuleRegistry.h" // JOBMODULE registration macro
#include "JobControl/JobCommand.h"
#include "MessageService/MsgService.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandEventListHandle.h"
#include "RecoBase/CandEventHandle.h"
#include "RecoBase/LinearFit.h"
#include "RecoBase/PropagationVelocity.h"
#include "CandTrackSR/CandTrackSRHandle.h"
#include "CandTrackSR/CandTrackSRListHandle.h"
#include "CandFitTrackSR/CandFitTrackSRHandle.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRFiducial.h"
#include "CandNtupleSR/NtpSRVertex.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CalDetSI/Helpers.h"
#include "CalDetDST/UberHit.h"
#include "CalDetDST/UberRecordLite.h"
#include "CalDetDST/UberRecHeader.h"
#include "CalDetDST/UberMC.h"
#include "CalDetDST/UberModuleLite.h"

#include "DatabaseInterface/DbiResultPtr.h"
#include "CalDetPID/CalDetBeamMomentum.h"
#include "Calibrator/CalTempCalibration.h"
#include "CalDetPID/CandCalDetPIDHandle.h"
#include "CalDetPID/NtpCalDetPID.h"

#include "Record/SimSnarlRecord.h"
#include "TParticle.h"

#include <cmath>

CVSID("$Id: UberModuleLite.cxx,v 1.19 2007/11/11 09:11:39 rhatcher Exp $");
JOBMODULE(UberModuleLite, "UberModuleLite", "A CalDet ntuple, based on NtupleBase");

using namespace std;

Int_t UberModuleLite::MakePlaneStripIndex(Int_t plane, Int_t strip)
{
   // ought to add range checking
   return plane*200 + strip;
   
}
// MAK: Feb 15, 2005,
// removed all vestiges of ancient fMC_MIP_Conversion, fMC_PE_Conversion hacks

//________________________________________________________________________________
UberModuleLite::UberModuleLite():
  fIsMC(kFALSE),
  fFirstEvent(kTRUE),
  fMCDecision(kFALSE),
  fRunNumber(0),
  fStarttime(0),
  candrec(0),
  fCDLH(0),
  fCDSI(0),
  fStripUidMap(),
  fShowerUidMap(),
  fTrackUidMap(),
  fEventUidMap()
{


  MSG("UberModuleLite",Msg::kDebug)<<"In UberModuleLite creator"<<endl;

}//end UberModuleLite()
//________________________________________________________________________________

UberModuleLite::~UberModuleLite()
{
  MSG("UberModuleLite",Msg::kDebug)<<"In ~UberModuleLite"<<endl;
}//end ~UberModuleLite()

//________________________________________________________________________________
JobCResult UberModuleLite::Get(MomNavigator* mom)
{
  MSG("UberModuleLite",Msg::kDebug)<<"In UberModuleLite get"<<endl;

  // Check that mom exists.
  assert(mom);  

  //ask mom  for a candidate record
  candrec=dynamic_cast<CandRecord*>
     (mom->GetFragment("CandRecord","PrimaryCandidateRecord"));
  if(candrec==0){
    MSG("UberModuleLite",Msg::kError)<<"candrec==0"<<endl;
    return JobCResult(JobCResult::kFailed);
  }
  if(fRunNumber==0){
    const CandHeader *ch = candrec->GetCandHeader();
    if(ch!=0){
      fRunNumber=ch->GetRun();
      MSG("UberModuleLite",Msg::kDebug)<<"RUN NUMBER "<<fRunNumber<<endl;
    }
  }

  fCDLH = dynamic_cast<CandDigitListHandle *>
    (candrec->FindCandHandle("CandDigitListHandle","canddigitlist"));
  if(fCDLH == 0) {//failure
    MSG("UberModuleLite",Msg::kError)<<"fCDLH==0"<<endl;
    return JobCResult(JobCResult::kFailed);
  }

  //now ask candrecord for a CandCalDetSIHandle
  fCDSI = dynamic_cast<CandCalDetSIHandle *>
    (candrec->FindCandHandle("CandCalDetSIHandle"));
  if(fCDSI == 0) {//failure
    MSG("UberModuleLite",Msg::kError)<<"fCDSI==0"<<endl;
    return JobCResult(JobCResult::kFailed);
  }

  
  MSG("UberModuleLite",Msg::kDebug)<<"Returning kPassed"<<endl;
  
  return JobCResult(JobCResult::kPassed); //success
}//end Get
//________________________________________________________________________________
JobCResult UberModuleLite::Reco(MomNavigator *mom)
{
  MSG("UberModuleLite",Msg::kDebug)<<"In UberModuleLite Reco"<<endl;

  //get an iterator over the digits
  CandDigitHandleItr hi(fCDLH->GetDaughterIterator());
  if(!hi.IsValid()){
    MSG("UberModuleLite",Msg::kError)<<"DaugherIter not valid"<<endl;
    return JobCResult(JobCResult::kFailed);
  }

  //Get a validity context and a plex handle
  const VldContext *vc = (*hi)->GetVldContext();
  PlexHandle ph(*vc, kTRUE);
  if(!fMCDecision){
    if(vc->GetSimFlag()==SimFlag::kReroot||vc->GetSimFlag()==SimFlag::kMC){
      fIsMC = kTRUE;
    }
    fMCDecision = kTRUE;
  }
  
  Calibrator::Instance().Reset(*vc);

  // Use validity context to create header for ntuple record
  UberRecHeader newHdr(*vc);
  newHdr.SetRunNo(fRunNumber);
  if(fFirstEvent){
    fStarttime = vc->GetTimeStamp().GetSec();
    
    fFirstEvent=kFALSE;
  }
  newHdr.SetStartTime(fStarttime);

  // Read temperature from the db.
  float temp = 100.0; // obviously fake value
  if(fIsMC){
       // set a bogus temperature
       newHdr.SetTemperature(temp);
  }
  else{
        newHdr.SetTemperature(Calibrator::Instance().GetTemperature());

  }

  // Read momentum from the db!
  float p = 0.1;// obviously fake value
  if(fIsMC){
       // try to figure out the momentum from mc truth
       // grab the simrecord
       const SimSnarlRecord* simrec=dynamic_cast<const SimSnarlRecord*>
            (mom->GetFragment("SimSnarlRecord"));
       if(simrec){
            // grab the stdhep block
            const TClonesArray* stdhep = dynamic_cast<const TClonesArray*> 
                 (simrec->FindComponent("TClonesArray", "StdHep"));
            if(stdhep){
                 // primary should be first in the list
                 const TParticle* part=
                      dynamic_cast<const TParticle*>(stdhep->At(0));
                 if(part) p = part->P();
            }
       }
  }
  else{
       DbiResultPtr<CalDetBeamMomentum> cdbm(*vc);
       if(cdbm.GetNumRows()==1) p = cdbm.GetRow(0)->GetMomentum();
  }
  newHdr.SetBeamMomentum(p);

  //create a new uberrecord
  UberRecordLite *snarldata = new UberRecordLite(newHdr);


//  MSG("UberModuleLite",Msg::kDebug)<<"Resetting Event"<<endl;
//  snarldata->ResetEvent();
  fStripUidMap.clear();
  fShowerUidMap.clear();
  fTrackUidMap.clear();
  fEventUidMap.clear();

  //start filling the snarldata

  //get the time since beginning of timeframe from from candigitlisthandle
  snarldata->triggertime = fCDLH->GetAbsTime();

  //get pid info from the CandCalDetSIHandle
  snarldata->snarlno = fCDSI->GetSnarl();
  MSG("UberModuleLite",Msg::kDebug)<<"SNARL NO "<<snarldata->snarlno<<endl;
  snarldata->triggerword = fCDSI->GetTrigSource();
  snarldata->ceradc[0]=fCDSI->GetKovADC2();
  snarldata->ceradc[1]=fCDSI->GetKovADC1();
  snarldata->ceradc[2]=fCDSI->GetKovADC3();
  snarldata->torbits=fCDSI->GetTriggerORBits();
  snarldata->torok=fCDSI->GetTriggerOROK();

  //get time relative to trigger time of cerenkov hits
  //note:  CalDetConstants::FARTIMECONVERT()--magic number to convert ticks to 
  //nanoseconds for the FAR DETECTOR electronics
  //we will need to fix this if the cerenkov is ever read out
  //using Near detector electronics (ie 2003)
  if(fCDSI->GetKovTimeStamp2()!=0){
    snarldata->certime[0]=fCDSI->GetKovTimeStamp2()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
    MSG("UberModuleLite",Msg::kDebug)<<"Kov 2 timestamp "<<fCDSI->GetKovTimeStamp2()
                              <<" time convert "<<CalDetConstants::FARTIMECONVERT
                              <<" trigger time "<<snarldata->triggertime
                              <<" in ns "
                              <<snarldata->triggertime/Munits::nanosecond
                              <<endl;
  }
  if(fCDSI->GetKovTimeStamp1()!=0){
    snarldata->certime[1]=fCDSI->GetKovTimeStamp1()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
  }
  if(fCDSI->GetKovTimeStamp3()!=0){
    snarldata->certime[2]=fCDSI->GetKovTimeStamp3()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
  }
  snarldata->toftdc[0]=fCDSI->GetTofTDC0();
  snarldata->toftdc[1]=fCDSI->GetTofTDC1();
  snarldata->toftdc[2]=fCDSI->GetTofTDC2();
  snarldata->tofadc[0]=fCDSI->GetTofADC0();
  snarldata->tofadc[1]=fCDSI->GetTofADC1();
  snarldata->tofadc[2]=fCDSI->GetTofADC2();
  if(fCDSI->GetTofTimeStamp()!=0){
     snarldata->toftime=fCDSI->GetTofTimeStamp()
        *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
     MSG("UberModuleLite",Msg::kDebug)<<"tof timestamp "<<fCDSI->GetTofTimeStamp()
                                      <<" time convert "<<CalDetConstants::FARTIMECONVERT
                                      <<setprecision(10)<<" trigger time "
                                      <<snarldata->triggertime
                                      <<" in ns "
                                      <<setprecision(10)
                                      <<snarldata->triggertime/Munits::nanosecond
                                      <<" tof time "<<snarldata->toftime<<endl;
     
  }
  if(fCDSI->GetTofADCTimeStamp0()!=0){
    snarldata->tofhittime[0]=fCDSI->GetTofADCTimeStamp0()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
  }
  if(fCDSI->GetTofADCTimeStamp1()!=0){
    snarldata->tofhittime[1]=fCDSI->GetTofADCTimeStamp1()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
  }
  if(fCDSI->GetTofADCTimeStamp2()!=0){
    snarldata->tofhittime[2]=fCDSI->GetTofADCTimeStamp2()
      *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
  }

  //add the number of seconds since the beginning of the run to the
  //trigger time to make trigger time absolute since beginning of run
  snarldata->triggertime+=vc->GetTimeStamp().GetSec()-fStarttime;
  MSG("UberModuleLite",Msg::kDebug)<<"Trigger time "<<snarldata->triggertime<<endl
                            <<"vld context time "
                            <<vc->GetTimeStamp().GetSec()<<endl;
/*
  PlexHandle ph(*vc, kTRUE);
  if(!fMCDecision){
    if(vc->GetSimFlag()==SimFlag::kReroot||vc->GetSimFlag()==SimFlag::kMC){
      fIsMC = kTRUE;
    }
    fMCDecision = kTRUE;
  }
*/
  Bool_t gotsigcorrconv = kFALSE;


  //check dead chip list in the CandCalDetSIHandle to find mindeadplane
  //and and number of dead planes
  set<UShort_t> deadplanemap;
  vector<RawChannelId>::const_iterator dcit(fCDSI->GetDeadChips().begin());
  MSG("UberModuleLite", Msg::kDebug)<<"Starting to get dead chip vector"<<endl
                               <<"Size of Dead Chip vector "
                               <<fCDSI->GetDeadChips().size()<<endl;
  while(dcit!=fCDSI->GetDeadChips().end()){
    //make a rawchannel id out of the dead chip (use va channel 3)
//changed 8-27-03 to make dead chip checking work with QIE chips
    RawChannelId newid;
    newid.SetDetector((*dcit).GetDetector());
    newid.SetElecType((*dcit).GetElecType()); 
    newid.SetCrate((*dcit).GetCrate());
    if((*dcit).GetElecType()==ElecType::kVA){
       newid.SetVarcId((*dcit).GetVarcId());
       newid.SetVmm((*dcit).GetVmm());
       newid.SetVaAdcSel((*dcit).GetVaAdcSel());
       newid.SetVaChip((*dcit).GetVaChip());
       newid.SetVaChannel(3);
    }
    if((*dcit).GetElecType()==ElecType::kQIE){
       newid.SetGeographicAddress((*dcit).GetGeographicAddress());
       newid.SetMasterChannel((*dcit).GetMasterChannel());
       newid.SetMinderChannel((*dcit).GetMinderChannel());
    }
//    newid.SetModeBits(false,false,false);
    newid.SetPedMode(false);
    newid.SetSparsMode(false);
    newid.SetCommonMode(false);
       

       
    if(ph.GetSEIdAltL(newid).GetPlane()>=0){
       deadplanemap.insert(ph.GetSEIdAltL(newid).GetPlane());
    }
    if(newid.GetElecType()==ElecType::kVA){
       if(fCDSI->GetCerenkovChannel3().GetElecType()==ElecType::kVA){
          if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel3())){
             deadplanemap.insert(65);
          }
       }
       if(fCDSI->GetCerenkovChannel1().GetElecType()==ElecType::kVA){
          if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel1())){
             deadplanemap.insert(66);
          }
       }
       if(fCDSI->GetCerenkovChannel2().GetElecType()==ElecType::kVA){
          if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel2())){
             deadplanemap.insert(67);
          }
       }
    }
    else if(newid.GetElecType()==ElecType::kQIE){
       if(fCDSI->GetCerenkovChannel3().GetElecType()==ElecType::kQIE){
          if(newid.IsSameChannel(fCDSI->GetCerenkovChannel3())){
             deadplanemap.insert(65);
          }
       }
       if(fCDSI->GetCerenkovChannel1().GetElecType()==ElecType::kQIE){
          if(newid.IsSameChannel(fCDSI->GetCerenkovChannel1())){
             deadplanemap.insert(66);
          }
       }
       if(fCDSI->GetCerenkovChannel2().GetElecType()==ElecType::kQIE){
          if(newid.IsSameChannel(fCDSI->GetCerenkovChannel2())){
             deadplanemap.insert(67);
          }
       }
    }
    //make another possible rawchannel id out of the dead chip(use va channel 11)
    //which could possibly be on the next plane
    //modified 8-17-03 to work with qie elec., 
    //need to look at 4 different masterchans to make sure we get all the planes
    RawChannelId newid2;
    newid2.SetDetector((*dcit).GetDetector());
    newid2.SetElecType((*dcit).GetElecType()); 
    newid2.SetCrate((*dcit).GetCrate());
    if((*dcit).GetElecType()==ElecType::kVA){
       newid2.SetVarcId((*dcit).GetVarcId());
       newid2.SetVmm((*dcit).GetVmm());
       newid2.SetVaAdcSel((*dcit).GetVaAdcSel());
       newid2.SetVaChip((*dcit).GetVaChip());
       newid2.SetVaChannel(11);
    }
    if((*dcit).GetElecType()==ElecType::kQIE){
       newid2.SetGeographicAddress((*dcit).GetGeographicAddress());
       newid2.SetMasterChannel((*dcit).GetMasterChannel()+1);
       newid2.SetMinderChannel((*dcit).GetMinderChannel());
    }
//    newid2.SetModeBits(false,false,false);
      newid2.SetPedMode(false);
      newid2.SetSparsMode(false);
      newid2.SetCommonMode(false);
       
  
    if(ph.GetSEIdAltL(newid2).GetPlane()>=0){
      deadplanemap.insert(ph.GetSEIdAltL(newid2).GetPlane());
    }
    if((*dcit).GetElecType()==ElecType::kQIE){
       RawChannelId newid3;
       newid3.SetDetector((*dcit).GetDetector());
       newid3.SetElecType((*dcit).GetElecType()); 
       newid3.SetCrate((*dcit).GetCrate());
       newid3.SetElecType((*dcit).GetElecType()); 
       newid3.SetGeographicAddress((*dcit).GetGeographicAddress());
       newid3.SetMasterChannel((*dcit).GetMasterChannel()+2);
       newid3.SetMinderChannel((*dcit).GetMinderChannel());
//       newid3.SetModeBits(false,false,false);
       newid3.SetPedMode(false);
       newid3.SetSparsMode(false);
       newid3.SetCommonMode(false);
       

       if(ph.GetSEIdAltL(newid3).GetPlane()>=0){
          deadplanemap.insert(ph.GetSEIdAltL(newid3).GetPlane());
       }
       RawChannelId newid4;
       newid4.SetDetector((*dcit).GetDetector());
       newid4.SetElecType((*dcit).GetElecType()); 
       newid4.SetCrate((*dcit).GetCrate());
       newid4.SetElecType((*dcit).GetElecType()); 
       newid4.SetGeographicAddress((*dcit).GetGeographicAddress());
       newid4.SetMasterChannel((*dcit).GetMasterChannel()+3);
       newid4.SetMinderChannel((*dcit).GetMinderChannel());
//       newid4.SetModeBits(false,false,false);
       newid4.SetPedMode(false);
       newid4.SetSparsMode(false);
       newid4.SetCommonMode(false);

    
       if(ph.GetSEIdAltL(newid4).GetPlane()>=0){
          deadplanemap.insert(ph.GetSEIdAltL(newid4).GetPlane());
       }
    }


    dcit++;
  }
  if(deadplanemap.size()!=0){
    snarldata->mindeadplaneno = *(deadplanemap.begin());
    snarldata->ndeadplanes = deadplanemap.size();
  }
  MSG("UberModuleLite",Msg::kDebug)<<"First dead plane "<<snarldata->mindeadplaneno
                            <<" No. dead planes "<<snarldata->ndeadplanes<<endl;
  


  //sort hits in event by plane
  CandDigitHandleKeyFunc* kf = hi.CreateKeyFunc();
  kf->SetFun(CalHelpers::KeyFromPlane);
  hi.GetSet()->AdoptSortKeyFunc(kf,kTRUE,kFALSE);
  //further sort hits by strip
  kf=hi.CreateKeyFunc();
  kf->SetFun(CalHelpers::KeyFromStrip);
  hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kFALSE);
  //and finally, sort by strip end
  kf=hi.CreateKeyFunc();
  kf->SetFun(CalHelpers::KeyFromEnd);
  hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kTRUE);
  kf = 0;
  MSG("UberModuleLite",Msg::kDebug)<<"Done sorting"<<endl;
  
  //declare variables to find maxes
  Float_t maxplanemip = 0.;
  Float_t maxp0stripmip=0.;
  Float_t maxp1stripmip=0.;
  Float_t maxp0hitmip=0.;
  Float_t maxp1hitmip=0.;
  Float_t totmipeven=0.;
  Float_t totmipodd=0.;

/*
  map<PlexStripEndId, REROOT_FLSDigit> FLSDigitMap;
  int nmchitstrips=0;
  if(fIsMC){
    const TClonesArray *hep = RerootExodus::GetStdHepList();
    mc->SetStdHep(hep);
    MSG("UberModuleLite",Msg::kDebug)<<"size of hep is "<<hep->GetEntries()<<endl;
    REROOT_Event *revt = gMINFast->GetREROOTEvent();
    mc->nhep = revt->n_stdheps();
    MSG("UberModuleLite",Msg::kDebug)<<"and mc nhep is  "<<mc->nhep<<endl;
    REROOT_StdHep *stdhep = static_cast<REROOT_StdHep *>(hep->At(0));
    if(stdhep!=0){
      mc->mctype = stdhep->ID();
      mc->mcpx = stdhep->Px();
      mc->mcpy = stdhep->Py();
      mc->mcpz = stdhep->Pz();
      mc->mcenergy = stdhep->E();
      mc->mcvx = stdhep->Xmm();
      mc->mcvy = stdhep->Ymm();
      mc->mcvz = stdhep->Zmm();
    }
    const TClonesArray *flsd = RerootExodus::GetFLSDigitList();
    int i=0;
    TObject *obj;
    while((obj=flsd->At(i))){
      REROOT_FLSDigit *mcdig = static_cast<REROOT_FLSDigit*>(obj);
      PlexStripEndId mcseid = RerootExodus::PECAB2SEId(mcdig->IPln(),
                                                       mcdig->IExtr(),
                                                       mcdig->ICell(),0);
      FLSDigitMap[mcseid]=*mcdig;
      i++;
    }
  }
*/
    

  
  //begin looping over planes
  MSG("UberModuleLite",Msg::kDebug)<<"begining to loop over hits"<<endl;;
  for(hi.Reset();hi.IsValid();hi.NextKey()){ 
    UShort_t plane = (*hi)->GetPlexSEIdAltL().GetPlane();
    MSG("UberModuleLite",Msg::kDebug)<<"PLANE "<<plane<<endl;
    Bool_t iscosmic = kFALSE;
    if(plane>=CalDetConstants::PLANECONST){ //if not a detector plane, continue
      iscosmic=kTRUE;
    }
    if(!iscosmic){
      snarldata->nhitplanes++; //increment nhitplanes
    }
    Float_t planemip=0.;

    //begin looping over strips
    for(CandDigitHandleItr siter(hi, kTRUE);
        siter.IsValid();siter.NextKey()){ 
      UShort_t strip = (*siter)->GetPlexSEIdAltL().GetBestSEId().GetStrip();
      MSG("UberModuleLite",Msg::kDebug)<<"STRIP "<<strip<<endl;
      if(!iscosmic){
        MSG("UberModuleLite",Msg::kDebug)<<"adding one to nhitstrips"<<endl;
        snarldata->nhitstrips++;  //increment nhitstrips
        MSG("UberModuleLite",Msg::kDebug)<<"adding next hit to tclonesarray"<<endl;
        //make a hit in UberEvent's TClonesArray
        int sindx = snarldata->AddNextHit(plane,strip);
        MSG("UberModuleLite",Msg::kDebug)<<"added next hit to tclonesarray"<<endl;
        //int uid = (*siter)->GetUidInt();
        Int_t uid = MakePlaneStripIndex(plane, strip);
        fStripUidMap.insert(std::make_pair(uid,sindx));

/*
        if(fIsMC){
          mc->AddNextHit(plane,strip);
        }
*/
      }
      if(iscosmic){
        snarldata->AddNextCosmicHit(plane,strip);
      }
      Float_t p0stripmip=0.;
      Float_t p1stripmip=0.;

      //begin looping over ends
      for(CandDigitHandleItr eiter(siter,kTRUE);
          eiter.IsValid();eiter.Next()){ //loop over ends
        
        MSG("UberModuleLite",Msg::kDebug)<<"RAW CHANNEL ID FOR THIS TIME ROUND:"<<endl;
//        (*eiter)->GetChannelId().Print();
        MSG("UberModuleLite",Msg::kDebug)<<endl;
//        (*eiter)->GetPlexSEIdAltL().Print();
        MSG("UberModuleLite",Msg::kDebug)<<" readouttype "
                                     <<ph.GetReadoutType((*eiter)->GetChannelId())
                                     <<endl;

        if(((*eiter)->GetPlexSEIdAltL()).size()!=1){
           MSG("UberModuleLite",Msg::kDebug)<<"Wrong size of Plexseidaltl "
                                        <<((*eiter)->GetPlexSEIdAltL()).size()
                                        <<" snarl no "<<snarldata->snarlno
                                        <<" plane is "<<plane<<endl;
        }
        if((*eiter)->GetPlexSEIdAltL().size()<=0){
          MSG("UberModuleLite",Msg::kDebug)<<"plexseidaltl is EMPTY"<<endl;
          //(*eiter)->GetChannelId().Print();
          //MSG("UberModuleLite",Msg::kDebug)<<endl;
          //(*eiter)->GetPlexSEIdAltL().Print();
          //MSG("UberModuleLite",Msg::kDebug)<<endl;
          continue;
        }
        /*
        //for debugging the reading of the near detector stuff:
        if((*eiter)->GetChannelId().GetElecType()==ElecType::kQIE){
          MSG("UberModuleLite",Msg::kDebug)<<" FOUND A QIE DIGIT "<<endl;
        }
        */
        PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
        StripEnd::StripEnd_t se = pse.GetEnd();
        PlexSEIdAltLItem pseitem = (*eiter)->GetPlexSEIdAltL().GetBestItem();


        //get hit info
        Int_t adc = (int)(*eiter)->GetCharge();
        Float_t siglin = (*eiter)->GetCharge(CalDigitType::kSigLin);
        Float_t npe = pseitem.GetPE();
        Float_t sigcorr = (*eiter)->GetCharge(CalDigitType::kSigCorr);
        Float_t mip = Calibrator::Instance().GetMIP((*eiter)
                                                   ->GetCharge(CalDigitType::kSigCorr),pse);
        Float_t time = ((*eiter)->GetSubtractedTime(CalTimeType::kT0)
                        /Munits::nanosecond);
        PlexPixelSpotId pps=ph.GetPixelSpotId(pse);
        Int_t agg = pps.GetEncoded();
        // added by mak March 25, 2003

        if(iscosmic){
          snarldata->AddNextCosmicHitValues(se, adc, time, agg);
          continue;
        }

        snarldata->nhits++;  //increment nhits
        if(!gotsigcorrconv){
          // MAK -- 6/3/03 : Added sig fpe protection
          if(mip>0) snarldata->sigcorrconv = sigcorr/mip;
          else snarldata->sigcorrconv=-1;
          gotsigcorrconv=kTRUE;
        }

        //add the hit info to UberEvent's TClonesArray
//CandDigit now removes the pedestal offset!!!!!, no need to do anything special
//for qie digits anymore
/*
        if((*eiter)->GetChannelId().GetElecType()==ElecType::kQIE){
          // MAK -- 6/3/03 : Added sig fpe protection
          float Nd=0;
          float Nc=0;
          if(siglin>0) Nd = 1.*adc/siglin;
          else Nd=1.0;
          if(mip) Nc = 1.*siglin/mip;
          else Nc=1.0;
          snarldata->AddNextHitValues(se, adc-50, siglin, 
                                      (adc-50)*1./(1.*adc/npe), (adc-50)/(Nd*Nc),
                                      time, agg);
        }
        else{
          snarldata->AddNextHitValues(se, adc, siglin, npe, mip, time, agg);
        }
*/

        snarldata->AddNextHitValues(se, adc, siglin, npe, mip, time, agg);
        MSG("UberModuleLite",Msg::kDebug)<<"STRIP END "<<se<<" adc "<<adc<<endl;
        planemip+=mip;  //increment the mips deposited in current plane

        if(pse.GetPlane()%2==0){  //calculate center in even planes
          snarldata->mipweighcentereven+=mip*pse.GetStrip();
          totmipeven+=mip;
        }
        else{ //calculate center in odd planes
          snarldata->mipweighcenterodd+=mip*pse.GetStrip();
          totmipodd+=mip;
        }

        //look for max hit adc, npe, mip, time
        if(adc>snarldata->maxadc){
          snarldata->maxadc = adc;
        }
        if(npe>snarldata->maxnpe){
          snarldata->maxnpe = npe;
        }
        if(mip>snarldata->maxmip){
          snarldata->maxmip = mip;
        }
        if(time>snarldata->maxtime){
          snarldata->maxtime = time;
        }

        //add mip in all planes but plane 0
        if(pse.GetPlane()!=0){
          snarldata->totmip+=mip;
        }

        //find mips in plane 0
        if(pse.GetPlane()==0){
          snarldata->p0totmip+=mip;
          p0stripmip+=mip;
          if(mip>maxp0hitmip){
            maxp0hitmip = mip;
            snarldata->p0maxmiptstamp = time;
          } 
        }

        //find mips in plane 1
        if(pse.GetPlane()==1){
          snarldata->p1totmip+=mip;
          p1stripmip+=mip;
          if(mip>maxp1hitmip){
            maxp1hitmip = mip;
            snarldata->p1maxmiptstamp = time;
          } 
        }
      }//end loop over ends
/*      
      if(fIsMC&&!iscosmic){
        PlexStripEndId mcpse = (*siter)->GetPlexSEIdAltL().GetBestSEId();
        mcpse.SetEnd(StripEnd::kNegative);
        REROOT_FLSDigit mcd = FLSDigitMap.find(mcpse)->second;
        MSG("UberModuleLite",Msg::kDebug)<<"pse.GetPlane() "
                                  <<(*siter)
          ->GetPlexSEIdAltL().GetBestSEId().GetPlane()
                                  <<" mc plane "<<mcd.IPln()<<endl;
        mc->AddNextHitValues(mcd.HitBits(),mcd.SumETrue(),mcd.TPos(),mcd.RawB(),
                             mcd.RawA(), mcd.CorrB(), mcd.CorrA(), mcd.TDCB(),
                             mcd.TDCA(), mcd.AveDistTrueB(),mcd.AveDistTrueA());
        mc->mcenergydep += mcd.SumETrue();
        nmchitstrips++;
      }
*/
      //find max mip in an strip in plane 0
      if(p0stripmip>maxp0stripmip){ 
        maxp0stripmip = p0stripmip;
        snarldata->p0stripmaxmip = strip;
      }
      //find max mip in a strip in plane 1
      if(p1stripmip>maxp1stripmip){
        maxp1stripmip = p1stripmip;
        snarldata->p1stripmaxmip = strip;
      }
    }//end loop over strips

    snarldata->mipweighaveplane+=planemip*plane;  //calculate mip weigh ave plane
    if(planemip>maxplanemip){  //find showermax
      maxplanemip = planemip;
      snarldata->mipshowermax = maxplanemip;
      snarldata->showermax = plane;
    }
  }//end loop over planes
  // MAK -- sept 21 2003
  if((snarldata->totmip+snarldata->p0totmip)>0){
       snarldata->mipweighaveplane/=(snarldata->totmip
                                     +snarldata->p0totmip);
  }
  else snarldata->mipweighaveplane=-1;
  // MAK -- 6/3/03
  // added some sig fpe checking here
  if(totmipeven>0) snarldata->mipweighcentereven/=totmipeven;
  else snarldata->mipweighcentereven=-1;
  if(totmipodd>0) snarldata->mipweighcenterodd/=totmipodd;
  else snarldata->mipweighcenterodd=-1;
  

  //loop again to calculate radius of event:
  for(hi.Reset();hi.IsValid();hi.NextKey()){
    UShort_t plane = ((*hi)->GetPlexSEIdAltL()).GetPlane();
    if(plane>=CalDetConstants::PLANECONST){
      continue;
    }
    float center = 0.;
    if(plane%2==0){
      center=snarldata->mipweighcentereven;
    }
    else{
      center = snarldata->mipweighcenterodd;
    }
    for(CandDigitHandleItr siter(hi, kTRUE);
        siter.IsValid();siter.NextKey()){ //loop over strips
      for(CandDigitHandleItr eiter(siter,kTRUE);
          eiter.IsValid();eiter.Next()){ //loop over ends
        PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
        if(pse.GetStrip()>=CalDetConstants::STRIPCONST){
          continue;
        }
        Float_t mip = Calibrator::Instance().GetMIP((*eiter)
                                                    ->GetCharge(CalDigitType::kSigCorr),pse);

        snarldata->mipweighrad += fabs(center-pse.GetStrip())*mip;
      }//end end loop
    }//end strip loop
  }//end planeloop

  // MAK -- 6/3/03 : Added sig fpe protection  
  if((snarldata->totmip+snarldata->p0totmip)>0.0) 
    snarldata->mipweighrad/=(snarldata->totmip+snarldata->p0totmip);
  else snarldata->mipweighrad=0;

  this->FillNtpCalDetPID(snarldata);

  this -> FillNtpShower(snarldata);
  this -> FillNtpTrack(snarldata);
  this -> FillNtpEvent(snarldata);

  MSG("UberModuleLite",Msg::kDebug)<<"giving fragment to mom"<<endl;  
  mom->AdoptFragment(snarldata);


  MSG("UberModuleLite",Msg::kDebug)<<"done with this event"<<endl;
  return JobCResult(JobCResult::kPassed);
}//end Reco()
//________________________________________________________________________________

void UberModuleLite::FillNtpCalDetPID(UberRecordLite* ntprec)
{
     MSG("UberModuleLite",Msg::kVerbose) 
          << "UberModuleLiteModule::FillNtpCalDetPID" << endl;

     const CandCalDetPIDHandle* pidh 
          = dynamic_cast <const CandCalDetPIDHandle*> 
          (candrec -> FindCandHandle("CandCalDetPIDHandle"));
     if ( !pidh ) return; // all done
     
     if(pidh->NoOverlap()) ntprec->cpid.nov=1;
     else ntprec->cpid.nov=0;

     if(pidh->InCERTime()) ntprec->cpid.inct=1;
     else ntprec->cpid.inct=0;
     
     ntprec->cpid.pid=pidh->GetPIDType();
     ntprec->cpid.olchi2=pidh->GetOLChi2();

     return;
}

void UberModuleLite::FillNtpShower(UberRecordLite* ntprec)
{
  //
  //  Purpose:  Private method used to fill shower portion of ntuple record.
  //
  //  Arguments: pointers to UberModuleLiteRecord and CandRecord
  //  
  //  Return: status.
  // 

  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLiteModule::FillNtpShower" << endl;

  const CandShowerListHandle *showerlisthandle 
   = dynamic_cast <const CandShowerListHandle*> 
     (candrec -> FindCandHandle("CandShowerListHandle"));
  if ( !showerlisthandle ) return; // all done

  Int_t nshower = 0;
  TIter showerItr(showerlisthandle->GetDaughterIterator());
  TClonesArray& showerarray = *(ntprec->shw);
  while (CandShowerHandle* shower=dynamic_cast<CandShowerHandle*>
                                                      (showerItr())) {
    Int_t uid = shower->GetUidInt();
    fShowerUidMap.insert(std::make_pair(uid,nshower));
    // Uses new with placement to efficiently create slice ntp
    NtpSRShower* ntpshower 
         = new(showerarray[nshower++])NtpSRShower(shower->GetNStrip());
    // Fill indices of associated strips in shower tree
    ntpshower->index = nshower - 1;
    ntpshower->ndigit = shower->GetNDigit();
    TIter showerstripItr(shower->GetDaughterIterator());
    Int_t nshowerstrip = 0;
//    while ( CandShowerHandle *showerstrip = dynamic_cast<CandShowerHandle*>
    //                                                (showerstripItr()) ) {
    while (TObject* obj = showerstripItr()) {
       CandStripHandle* showerstrip = dynamic_cast<CandStripHandle*>(obj);
       if(showerstrip){
          Int_t plane = showerstrip->GetPlane();
          Int_t strip = showerstrip->GetStrip();
//        Int_t uid = showerstrip->GetUidInt();
          Int_t uid = MakePlaneStripIndex(plane,strip);
          Int_t stripindex = fStripUidMap[uid]; // should use find
          ntpshower->AddStripAt(stripindex,nshowerstrip); // add index to strip
          nshowerstrip++;
       }
    }

    // Set range of planes included in slice
    ntpshower->plane.n = shower->GetNPlane();
    ntpshower->plane.nu = shower->GetNPlane(PlaneView::kU);
    ntpshower->plane.nv = shower->GetNPlane(PlaneView::kV);
    ntpshower->plane.beg = shower->GetBegPlane();
    ntpshower->plane.begu = shower->GetBegPlane(PlaneView::kU);
    ntpshower->plane.begv = shower->GetBegPlane(PlaneView::kV);
    ntpshower->plane.end = shower->GetEndPlane();
    ntpshower->plane.endu = shower->GetEndPlane(PlaneView::kU);
    ntpshower->plane.endv = shower->GetEndPlane(PlaneView::kV);
    // Set summed charge in shower
    ntpshower->ph.raw = shower->GetCharge(CalStripType::kNone);
    ntpshower->ph.siglin = shower->GetCharge(CalStripType::kSigLin);
    ntpshower->ph.sigcor = shower->GetCharge(CalStripType::kSigCorr);
    ntpshower->ph.pe = shower->GetCharge(CalStripType::kPE);
    ntpshower->ph.mip = shower->GetCharge(CalStripType::kMIP);
    ntpshower->ph.gev = shower->GetCharge(CalStripType::kGeV);

  }

  return;
}//end FillNtpShower(UberRecordLite, candrec)

//________________________________________________________________________________

void UberModuleLite::FillNtpTrack(UberRecordLite* ntprec)
{
  //
  //  Purpose:  Private method used to fill track portion of ntuple record.
  //
  //  Arguments: pointers to UberModuleLiteRecord and CandRecord
  //  
  //  Return: status.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrack" << endl;

  const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y

  const CandTrackListHandle *tracklisthandle 
     = dynamic_cast <const CandTrackListHandle*> 
     (candrec -> FindCandHandle("CandFitTrackListHandle"));
  if ( !tracklisthandle ) {
    tracklisthandle = dynamic_cast <const CandTrackListHandle*>
      (candrec -> FindCandHandle("CandTrackListHandle"));
  }
  if ( !tracklisthandle ) return;

  CandTrackSRListHandle* tracksrlisthandle=dynamic_cast<CandTrackSRListHandle*>
                      (candrec->FindCandHandle("CandTrackSRListHandle"));

  const VldContext& vld = *(candrec->GetVldContext());

  TIter trackItr(tracklisthandle->GetDaughterIterator());
  TClonesArray& trackarray = *(ntprec->trk);
  Int_t ntrack = 0;
  while (CandTrackHandle* track = dynamic_cast<CandTrackHandle*>(trackItr())){
    Int_t uid = track->GetUidInt();
    fTrackUidMap.insert(std::make_pair(uid,ntrack));
    // Uses new with placement to efficiently create event ntp
    NtpSRTrack* ntptrack 
             = new(trackarray[ntrack++])NtpSRTrack(track->GetNStrip());

    ntptrack->index = ntrack-1;

    CandTrackSRHandle *tracksr = dynamic_cast<CandTrackSRHandle*>(track);
    //CandFitTrackHandle *fittrack = dynamic_cast<CandFitTrackHandle*>(track);
    CandFitTrackSRHandle *fittracksr
                                =dynamic_cast<CandFitTrackSRHandle*>(track);

    // Set range of planes included in track
    NtpSRTrackPlane& plane = ntptrack->plane;
    plane.n = track->GetNPlane();
    plane.nu = track->GetNPlane(PlaneView::kU);
    plane.nv = track->GetNPlane(PlaneView::kV);
    plane.beg = track->GetBegPlane();
    plane.begu = track->GetBegPlane(PlaneView::kU);
    plane.begv = track->GetBegPlane(PlaneView::kV);
    plane.end = track->GetEndPlane();
    plane.endu = track->GetEndPlane(PlaneView::kU);
    plane.endv = track->GetEndPlane(PlaneView::kV);
    if (tracksr) plane.ntrklike = tracksr->GetNTrackPlane();
    if (fittracksr) plane.ntrklike = fittracksr->GetNTrackPlane();

    // Set summed pulse height information for track
    NtpSRStripPulseHeight& ph = ntptrack->ph;
    ph.raw = track->GetCharge(CalStripType::kNone);
    ph.siglin = track->GetCharge(CalStripType::kSigLin);
    ph.sigcor = track->GetCharge(CalStripType::kSigCorr);
    ph.pe = track->GetCharge(CalStripType::kPE);
    ph.mip = track->GetCharge(CalStripType::kMIP);
    ph.gev = track->GetCharge(CalStripType::kGeV);
    
    ntptrack->ds = track->GetdS(); // distance from vertex to end
    ntptrack->range = track->GetRange(); // g/cm**2 from vertex to end
    // CPU to create track list
    if (tracksrlisthandle) ntptrack->cputime = tracksrlisthandle->GetCPUTime();

    this->FillNtpTrackVertex(ntptrack,track);
    this->FillNtpTrackEnd(ntptrack,track);
    this->FillNtpTrackLinearFit(ntptrack,track);

    this->FillNtpTrackFidVtx(ntptrack,track,vld);
    this->FillNtpTrackFidEnd(ntptrack,track,vld);
    this->FillNtpTrackFidAll(ntptrack,track,vld);

    this->FillNtpTrackMomentum(ntptrack,track);
    this->FillNtpTrackFit(ntptrack,track);
    this->FillNtpTrackTime(ntptrack,track);

    // Loop over strips associated with track
    TIter trackstripItr(track->GetDaughterIterator());
    Int_t ntrackstrip = 0;
    // Fill indices of associated strips in track ntuple
//    while ( CandStripHandle *trackstrip = dynamic_cast<CandStripHandle*>
    //                                              (trackstripItr())) {
    while (TObject* obj = (trackstripItr())) {
       CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>(obj);
       if(trackstrip){
          Int_t plane = trackstrip->GetPlane();
          Int_t strip = trackstrip->GetStrip();
          Int_t uid = MakePlaneStripIndex(plane,strip);
          Int_t stripindex = fStripUidMap[uid];
          ntptrack->AddStripAt(stripindex,ntrackstrip); // add index to track
          // Track positional information at plane associated with each strip
          Int_t iplane = trackstrip->GetPlane();
          ntptrack->stpu[ntrackstrip] = track->GetU(iplane);
          ntptrack->stpv[ntrackstrip] = track->GetV(iplane);
          ntptrack->stpx[ntrackstrip] = cos45*(ntptrack->stpu[ntrackstrip]
                                               -ntptrack->stpv[ntrackstrip]);
          ntptrack->stpy[ntrackstrip] = cos45*(ntptrack->stpu[ntrackstrip]
                                               +ntptrack->stpv[ntrackstrip]);
          ntptrack->stpz[ntrackstrip] = track->GetZ(iplane);
          // dS is travel distance from vertex
          ntptrack->stpds[ntrackstrip] = track->GetdS(iplane);
          
          // Fit track dependent quantities
          if ( fittracksr ) {
             ntptrack->stpfitchi2[ntrackstrip] = fittracksr->GetPlaneChi2(iplane);
             ntptrack->stpfitprechi2[ntrackstrip] 
                = fittracksr->GetPlanePreChi2(iplane);
             ntptrack->stpfitqp[ntrackstrip] = fittracksr->GetPlaneQP(iplane);
             if ( fittracksr->GetBadFit(iplane) ) ntptrack->stpfit[ntrackstrip] = 0;
          }

          // Strip end dependent quantities
          for (UInt_t i = 0; i < 2; i++ ) {
             StripEnd::EStripEnd stripend = StripEnd::kNegative;
             if (i == 1) stripend = StripEnd::kPositive;
             if ( trackstrip->GetNDigit(stripend) > 0 ) {
                Float_t sigmap = track->GetStripCharge(trackstrip,
                                                       CalStripType::kSigMapped,stripend);
                Float_t mip = track->GetStripCharge(trackstrip,
                                                    CalStripType::kMIP,stripend);
                Float_t gev = track->GetStripCharge(trackstrip,
                                                    CalStripType::kGeV,stripend);
                ntptrack->SetPh(sigmap,mip,gev,ntrackstrip,i); 
                ntptrack->SetTime(track->GetT(iplane,stripend),ntrackstrip,i);
             }
             // NJT 07/04
             // Use the new Calibrator to dig out the values of interest.
             float c0 = Calibrator::Instance().DecalStripToStrip(1.0,trackstrip->GetStripEndId(stripend));
             ntptrack->SetAttnC0( c0 ,ntrackstrip,i);

             float t0 = Calibrator::Instance().DecalTime(0.0,trackstrip->GetCharge(CalDigitType::kNone,stripend),trackstrip->GetStripEndId(stripend));
             ntptrack->SetCalT0(t0,ntrackstrip,i);

          }
          ntrackstrip++;
       }
    }
        
  }

  return;

}//end FillNtpTrack()

//________________________________________________________________________________

void UberModuleLite::FillNtpEvent(UberRecordLite* ntprec)
{
  //
  //  Purpose:  Private method used to fill event portion of ntuple record.
  //
  //  Arguments: pointers to UberModuleLiteRecord and CandRecord
  //  
  //  Return: status.
  // 

  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpEvent" << endl;

  const CandEventListHandle *eventlisthandle 
   = dynamic_cast <const CandEventListHandle*> 
     (candrec -> FindCandHandle("CandEventListHandle"));
  if ( !eventlisthandle ) return; // all done

  Int_t nevent = 0;
  TIter eventItr(eventlisthandle->GetDaughterIterator());
  TClonesArray& eventarray = *(ntprec->evt);
  while (CandEventHandle* event=dynamic_cast<CandEventHandle*>
                                                      (eventItr())) {
    Int_t uid = event->GetUidInt();
    fEventUidMap.insert(std::make_pair(uid,nevent));
    // Uses new with placement to efficiently create event ntp
    NtpSREvent* ntpevent 
    = new(eventarray[nevent++])NtpSREvent(event->GetNStrip(),
                      event->GetLastShower()+1,event->GetLastTrack()+1);
    ntpevent->index = nevent-1;
    ntpevent->ndigit = event->GetNDigit();

    TIter eventstripItr(event->GetDaughterIterator());
    Int_t neventstrip = 0;
    // Fill indices of associated strips,showers,tracks in event ntuple
    while ( CandStripHandle *eventstrip = dynamic_cast<CandStripHandle*>
                                                      (eventstripItr())) {
      Int_t stripindex = fStripUidMap[eventstrip->GetUidInt()];
      ntpevent->AddStripAt(stripindex,neventstrip); // add index to event
      neventstrip++;
    }
    for (Int_t i = 0; i <= event->GetLastTrack(); i++ ) {
      const CandTrackHandle* track = event->GetTrack(i);
      Int_t trackindex = fTrackUidMap[track->GetUidInt()];
      ntpevent -> AddTrackAt(trackindex,i);
    }
    for (Int_t i = 0; i <= event->GetLastShower(); i++ ) {
      const CandShowerHandle* shower = event->GetShower(i);
      Int_t showerindex = fShowerUidMap[shower->GetUidInt()];
      ntpevent -> AddShowerAt(showerindex,i);
    }

    // Set range of planes included in event
    ntpevent->plane.n = event->GetNPlane();
    ntpevent->plane.nu = event->GetNPlane(PlaneView::kU);
    ntpevent->plane.nv = event->GetNPlane(PlaneView::kV);
    ntpevent->plane.beg = event->GetBegPlane();
    ntpevent->plane.begu = event->GetBegPlane(PlaneView::kU);
    ntpevent->plane.begv = event->GetBegPlane(PlaneView::kV);
    ntpevent->plane.end = event->GetEndPlane();
    ntpevent->plane.endu = event->GetEndPlane(PlaneView::kU);
    ntpevent->plane.endv = event->GetEndPlane(PlaneView::kV);
    // Set summed charge in event
    ntpevent->ph.raw = event->GetCharge(CalStripType::kNone);
    ntpevent->ph.siglin = event->GetCharge(CalStripType::kSigLin);
    ntpevent->ph.sigcor = event->GetCharge(CalStripType::kSigCorr);
    ntpevent->ph.pe = event->GetCharge(CalStripType::kPE);
    ntpevent->ph.mip = event->GetCharge(CalStripType::kMIP);
    ntpevent->ph.gev = event->GetCharge(CalStripType::kGeV);
  }

  return;
}//end FillNtpEvent()
//________________________________________________________________________________


void UberModuleLite::FillNtpTrackMomentum(NtpSRTrack* ntptrack,
                                       const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack momentum data member
  //            given a pointer to the ntuple track and a pointer to the
  //            associated CandTrackHandle.
  //
  //  Return: none.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackMomentum" << endl;

  NtpSRMomentum& momentum = ntptrack->momentum;
  momentum.range = track->GetMomentum();

  const CandFitTrackHandle* fittrack 
    = dynamic_cast<const CandFitTrackHandle*>(track);
  if ( fittrack ) {
      if(fittrack->GetMomentumCurve())
          momentum.qp = fittrack->GetEMCharge()/fittrack->GetMomentumCurve();
      else momentum.qp=0;
  }

  const CandFitTrackSRHandle* fittracksr 
    = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {
    momentum.eqp = fittracksr->GetVtxQPError();
  }

  return;

}

void UberModuleLite::FillNtpTrackTime(NtpSRTrack* ntptrack,
                                   const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack time data member
  //            given a pointer to the ntuple track and a pointer to the
  //            associated CandTrackHandle.
  //
  //  Return: none.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackTime" << endl;

  NtpSRTrackTime& time = ntptrack->time;

  const CandTrackSRHandle* tracksr 
    = dynamic_cast<const CandTrackSRHandle*>(track);
  if ( tracksr ) {
    time.ndigit = tracksr->GetNTimeFitDigit();
    time.chi2   = tracksr->GetTimeFitChi2();
  }

  const CandFitTrackSRHandle* fittracksr 
    = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {
    time.ndigit = fittracksr->GetNTimeFitDigit();
    time.chi2   = fittracksr->GetTimeFitChi2();
  }

  time.cdtds = fabs(track->GetTimeSlope())*3.e8;
  time.dtds = track->GetTimeSlope();
  time.t0   = track->GetTimeOffset();
  
  Double_t totuph[2] = {0}, totvph[2] = {0};
  Int_t    ndutimespace[1000] = {0}, ndvtimespace[1000] = {0};
  Double_t dutimespace[1000]  = {0}, dvtimespace[1000] = {0};
  
  TIter trackstripItr(track->GetDaughterIterator());
  while (CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
                                                   (trackstripItr())) {
    Int_t iplane = trackstrip -> GetPlane();
    assert(iplane >= 0 && iplane < 1000);
    Float_t ph0 = trackstrip -> GetCharge(StripEnd::kNegative);
    Float_t ph1 = trackstrip -> GetCharge(StripEnd::kPositive);
    Float_t t0  = track->GetT(iplane,StripEnd::kNegative);
    Float_t t1  = track->GetT(iplane,StripEnd::kPositive);

    if ( trackstrip -> GetPlaneView() == PlaneView::kU ) {
      if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
        totuph[0] += ph0;
        time.u0   += ph0*t0;
      }
      if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
        totuph[1] += ph1;
        time.u1   += ph1*t1;
      }
      if ( t0 > 0 && t1 > 0 ) {
        Float_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
        if (!ndutimespace[iplane] || fabs(dapos)<fabs(dutimespace[iplane])) {
          dutimespace[iplane] = dapos;
          ndutimespace[iplane] = 1;
        }
      }
    }

    if ( trackstrip -> GetPlaneView() == PlaneView::kV ) {
      if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
        totvph[0] += ph0;
        time.v0   += ph0*t0;
      }
      if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
        totvph[1] += ph1;
        time.v1   += ph1*t1;
      }
      if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 &&
           trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
        Float_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
        if (!ndvtimespace[iplane]||fabs(dapos) < fabs(dvtimespace[iplane])) {
          dvtimespace[iplane] = dapos;
          ndvtimespace[iplane] = 1;
        }
      }
    }
  }

  Int_t ndu = 0;
  Int_t ndv = 0;
  for (int ip = 0; ip < 1000; ip++ ) {
    if ( ndutimespace[ip] > 0 ) {
      ndu++;
      time.du += dutimespace[ip];
    }
    if ( ndvtimespace[ip] > 0 ) {
      ndv++;
      time.dv += dvtimespace[ip];
    }
  }
  if ( totuph[0] > 0. ) time.u0 /= totuph[0];
  if ( totvph[0] > 0. ) time.v0 /= totvph[0];
  if ( totuph[1] > 0. ) time.u1 /= totuph[1];
  if ( totvph[1] > 0. ) time.v1 /= totvph[1];
  if ( ndu ) time.du /= (Float_t)ndu;
  if ( ndv ) time.dv /= (Float_t)ndv;

  return;

}

void UberModuleLite::FillNtpTrackFit(NtpSRTrack* ntptrack,
                                  const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack fit data member
  //            given a pointer to the ntuple track and a pointer to the
  //            associated CandTrackHandle.
  //
  //  Return: none.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackFit" << endl;

  NtpSRFitTrack& fit = ntptrack->fit;

  const CandFitTrackHandle* fittrack 
    = dynamic_cast<const CandFitTrackHandle*>(track);
  if ( fittrack ) {
    fit.chi2 = fittrack->GetChi2();
    fit.pass = fittrack->GetPass();
  }

  const CandFitTrackSRHandle* fittracksr 
    = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {
    fit.ndof      = fittracksr->GetNDOF();
    fit.niterate  = fittracksr->GetNIterate();
    fit.nswimfail = fittracksr->GetNSwimFail();
    fit.cputime   = fittracksr->GetCPUTime();
  }

  return;

}

void UberModuleLite::FillNtpTrackFidVtx(NtpSRTrack* ntptrack,
                                     const CandTrackHandle* track,
                                     const VldContext& vld) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack fidvtx data member
  //            given a pointer to the ntuple track, a pointer to the
  //            associated CandTrackHandle, and the vldcontext of the
  //            event.
  //
  //  Return: none.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackFidVtx" << endl;

  NtpSRFiducial& fid = ntptrack->fidvtx;
  NtpSRVertex& vtx = ntptrack->vtx;
  //NtpSRFitVertex& vtx = ntptrack->vtx;
  this->FillNtpFiducial(fid,vtx,vld); // fills dr,dz

  fid.trace   = track->GetVtxTrace();
  fid.tracez  = track->GetVtxTraceZ();

  const CandFitTrackSRHandle* fittracksr 
    = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {                                     
    fid.nplane  = fittracksr->GetVtxExtrapolate();
    //   fid.nplaneu = fittracksr->GetVtxExtrapolate(PlaneView::kU);
    //  fid.nplanev = fittracksr->GetVtxExtrapolate(PlaneView::kV);
  }  

  return;

}

void UberModuleLite::FillNtpTrackFidEnd(NtpSRTrack* ntptrack,
                                     const CandTrackHandle* track,
                                     const VldContext& vld) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack fidend data member
  //            given a pointer to the ntuple track, a pointer to the
  //            associated CandTrackHandle, and the vldcontext of the
  //            event.
  //
  //  Return: none.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackFidEnd" << endl;

  NtpSRFiducial& fid = ntptrack->fidend;
  NtpSRVertex& end = ntptrack->end;
  //NtpSRFitVertex& end = ntptrack->end;
  this->FillNtpFiducial(fid,end,vld); // fills dr,dz

  fid.trace   = track->GetEndTrace();
  fid.tracez  = track->GetEndTraceZ();

  const CandFitTrackSRHandle* fittracksr 
    = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {                                     
    fid.nplane  = fittracksr->GetEndExtrapolate();
    //  fid.nplaneu = fittracksr->GetEndExtrapolate(PlaneView::kU);
    //  fid.nplanev = fittracksr->GetEndExtrapolate(PlaneView::kV);
  }  

  return;

}

void UberModuleLite::FillNtpTrackFidAll(NtpSRTrack* ntptrack,
                                     const CandTrackHandle* track,
                                     const VldContext& vld) {
  //
  //  Purpose:  Private method used to fill NtpSRTrack fidall data member
  //            given a pointer to the ntuple track, a pointer to the
  //            associated CandTrackHandle, and the vldcontext of the
  //            event.
  //
  //  Return: none.
  //
  //  Notes: This method should be called after FillNtpTrackFidAll
  //         and FillNtpFiducialEnd 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpFiducialAll" << endl;
  const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y

  NtpSRFiducial& fidall = ntptrack->fidall;

  NtpSRFiducial fidtmp;
  NtpSRVertex vtxtmp;
  fidall.dr = min(ntptrack->fidvtx.dr,ntptrack->fidend.dr);
  fidall.dz = min(ntptrack->fidvtx.dz,ntptrack->fidend.dz);
  fidall.trace   = min(ntptrack->fidvtx.trace,ntptrack->fidend.trace);
  fidall.tracez  = min(ntptrack->fidvtx.tracez,ntptrack->fidend.tracez);
  fidall.nplane  = min(ntptrack->fidvtx.nplane,ntptrack->fidend.nplane);
  //  fidall.nplaneu  = min(ntptrack->fidvtx.nplaneu,ntptrack->fidend.nplaneu);
  //  fidall.nplanev  = min(ntptrack->fidvtx.nplanev,ntptrack->fidend.nplanev);

  TIter trackstripItr(track->GetDaughterIterator());
  while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
        (trackstripItr())) {
    // For all strips calculate closest approach to boundaries
    Int_t iplane = trackstrip->GetPlane();
    vtxtmp.plane = iplane;
    vtxtmp.u = track->GetU(iplane);
    vtxtmp.v = track->GetV(iplane);
    vtxtmp.x = cos45*(vtxtmp.u-vtxtmp.v);
    vtxtmp.y = cos45*(vtxtmp.u+vtxtmp.v);
    this->FillNtpFiducial(fidtmp,vtxtmp,vld); // fills dr,dz
    fidall.dr = min(fidtmp.dr,fidall.dr);
    fidall.dz = min(fidtmp.dz,fidall.dz);
  }

  return;

}

void UberModuleLite::FillNtpTrackVertex(NtpSRTrack* ntptrack,
                                     const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill vertex portion of ntuple track.
  //
  //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
  //  
  //  Return: status.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackVertex" << endl;

  const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y

  NtpSRVertex& vtx = ntptrack->vtx;
  //  NtpSRFitVertex& vtx = ntptrack->vtx;

  vtx.u = track->GetVtxU();
  vtx.v = track->GetVtxV();
  vtx.x = cos45*(vtx.u-vtx.v);
  vtx.y = cos45*(vtx.u+vtx.v);
  vtx.z = track->GetVtxZ();
  vtx.t = track->GetVtxT();
  vtx.plane = track->GetVtxPlane();
  vtx.dcosu = track->GetDirCosU();
  vtx.dcosv = track->GetDirCosV();
  vtx.dcosx = cos45*(vtx.dcosu-vtx.dcosv);
  vtx.dcosy = cos45*(vtx.dcosu+vtx.dcosv);
  vtx.dcosz = track->GetDirCosZ();

  const CandFitTrackSRHandle* fittracksr
                        = dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {
    vtx.eu = fittracksr->GetVtxUError();
    vtx.ev = fittracksr->GetVtxVError();
    vtx.ex = cos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
    vtx.ey = cos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
    Double_t edudz = fittracksr->GetVtxdUError();
    Double_t edvdz = fittracksr->GetVtxdVError();
    // These calculations should include the dudz and dvdz covariance terms
    // but currently the covariance terms are not accessible
    vtx.edcosz = fabs(vtx.dcosz)*sqrt(pow(vtx.dcosu*edudz,2)+
                                      pow(vtx.dcosv*edvdz,2));
    vtx.edcosu = sqrt(fabs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edvdz,2)
               + pow((pow(vtx.dcosz,2)+pow(vtx.dcosv,2))*edudz,2));
    vtx.edcosv = sqrt(fabs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edudz,2)
               + pow((pow(vtx.dcosz,2)+pow(vtx.dcosu,2))*edvdz,2));
    vtx.edcosx = cos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
    vtx.edcosy = cos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
  }

  return;

}

void UberModuleLite::FillNtpTrackEnd(NtpSRTrack* ntptrack,
                                  const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill end portion of ntuple track.
  //
  //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
  //  
  //  Return: status.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpTrackEnd" << endl;

  const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y

  NtpSRVertex& end = ntptrack->end;
  //  NtpSRFitVertex& end = ntptrack->end;

  end.u = track->GetEndU();
  end.v = track->GetEndV();
  end.x = cos45*(end.u-end.v);
  end.y = cos45*(end.u+end.v);
  end.z = track->GetEndZ();
  end.t = track->GetEndT();
  end.plane = track->GetEndPlane();
  end.dcosu = track->GetEndDirCosU();
  end.dcosv = track->GetEndDirCosV();
  end.dcosx = cos45*(end.dcosu-end.dcosv);
  end.dcosy = cos45*(end.dcosu+end.dcosv);
  end.dcosz = track->GetEndDirCosZ();

  const CandFitTrackSRHandle* fittracksr
                    =dynamic_cast<const CandFitTrackSRHandle*>(track);
  if ( fittracksr ) {
    end.eu = fittracksr->GetEndUError();
    end.ev = fittracksr->GetEndVError();
    end.ex = cos45*sqrt(end.eu*end.eu+end.ev*end.ev);
    end.ey = cos45*sqrt(end.eu*end.eu+end.ev*end.ev);
    Double_t edudz = fittracksr->GetEnddUError();
    Double_t edvdz = fittracksr->GetEnddVError();
    // These calculations should include the dudz and dvdz covariance terms
    // but currently the covariance terms are not accessible
    end.edcosz = fabs(end.dcosz)*sqrt(pow(end.dcosu*edudz,2)+
                                      pow(end.dcosv*edvdz,2));
    end.edcosu = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edvdz,2)
               + pow((pow(end.dcosz,2)+pow(end.dcosv,2))*edudz,2));
    end.edcosv = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edudz,2)
               + pow((pow(end.dcosz,2)+pow(end.dcosu,2))*edvdz,2));
    end.edcosx = cos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
    end.edcosy = cos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
  }

  return;

}

void UberModuleLite::FillNtpTrackLinearFit(NtpSRTrack* ntptrack,
                                        const CandTrackHandle* track) {
  //
  //  Purpose:  Private method used to fill linar fit track portion of 
  //            ntuple track.
  //
  //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
  //  
  //  Return: status.
  // 


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpFitTrack" << endl;

  const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y

  // First array dimension is view (u or v)
  // Second array dimension is plane number
  const Int_t nplane = 1000;
  const Int_t nview = 2; // u,v
  Double_t zfit[nview][nplane]={{0},{0}},fit[nview][nplane]={{0},{0}};
  Double_t wfit[nview][nplane]={{0},{0}},ph[nview][nplane]={{0},{0}};

  // Loop over all strips on track
  TIter trackstripItr(track->GetDaughterIterator());
  while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
        (trackstripItr())) {

    Float_t tpos = trackstrip->GetTPos();
    Float_t phend[2]  = {trackstrip->GetCharge(StripEnd::kNegative),
                         trackstrip->GetCharge(StripEnd::kPositive)};
    Float_t phsum   = phend[0] + phend[1];

    Int_t iplane = trackstrip->GetPlane();
    assert(iplane >= 0 && iplane < 1000);

    Int_t iview = -1;
    if ( trackstrip->GetPlaneView()==PlaneView::kU ) iview = 0;
    else if ( trackstrip->GetPlaneView()==PlaneView::kV ) iview = 1;
    else continue;

    zfit[iview][iplane] = trackstrip->GetZPos();
    fit[iview][iplane] += tpos*phsum; // pulse height weighted average of tpos
    ph[iview][iplane]  += phsum;  // summed pulse height both ends 
    wfit[iview][iplane] = 1;   // weight assigned for use in fit
  }

  // Finished with loop over all strips, calculate fit results
  for ( Int_t ip = 0; ip < nplane; ip++ ) {
    for ( Int_t iv = 0; iv < nview; iv++ ) {
      if ( ph[iv][ip] > 0 ) fit[iv][ip] /= ph[iv][ip];
    }
  }

  Double_t uparm[2],ueparm[2],vparm[2],veparm[2];
  LinearFit::Weighted(1000,zfit[0],fit[0],wfit[0],uparm,ueparm);
  LinearFit::Weighted(1000,zfit[1],fit[1],wfit[1],vparm,veparm);
  NtpSRVertex& lin = ntptrack->lin;

  Double_t dudz = uparm[1];
  Double_t dvdz = vparm[1];
  Double_t dxdz = cos45*(uparm[1]-vparm[1]);
  Double_t dydz = cos45*(uparm[1]+vparm[1]);
  lin.z = track->GetVtxZ();
  lin.t = track->GetVtxT();
  lin.u = uparm[0]+lin.z*dudz;
  lin.v = vparm[0]+lin.z*dvdz;
  lin.x = cos45*(uparm[0]-vparm[0])+lin.z*dxdz;
  lin.y = cos45*(uparm[0]+vparm[0])+lin.z*dydz;
  Double_t anorm = sqrt(1.+dxdz*dxdz+dydz*dydz);
  Int_t sign = 1;
  if ( track->GetTimeSlope() < 0. ) sign = -1;
  lin.dcosu = sign*dudz/anorm;
  lin.dcosv = sign*dvdz/anorm;
  lin.dcosx = sign*dxdz/anorm;
  lin.dcosy = sign*dydz/anorm;
  lin.dcosz = sign*1./anorm;

  return;

}





void UberModuleLite::FillNtpFiducial(NtpSRFiducial& fid, const NtpSRVertex& vtx, 
                                  const VldContext& vld) {
  //
  //  Purpose:  Private method used to fill NtpSRFiducial given a point
  //            defined by NtpSRVertex and the detector type.
  //
  //  Return: none.
  //  


  MSG("UberModuleLite",Msg::kVerbose) << "UberModuleLite::FillNtpFiducial" << endl;

  UgliGeomHandle ugh(vld);
  Float_t zextent[2];
  ugh.GetZExtent(zextent[0],zextent[1]); // apparently this doesn't work yet?
  Detector::Detector_t dtype = vld.GetDetector();

  Float_t du,dv,dx,dy;

  switch (dtype) {

  case Detector::kFar:
    du = min(4.-vtx.u,4.+vtx.u);
    dv = min(4.-vtx.v,4.+vtx.v);
    dx = min(4.-vtx.x,4.+vtx.x);
    dy = min(4.-vtx.y,4.+vtx.y);
    fid.dr = min(min(du,dv),min(dx,dy));
    if ( fid.dr < 0. ) fid.dr = 0;
    fid.dz = min(vtx.z-zextent[0],zextent[1]-vtx.z);
    
    break;

  case Detector::kCalDet:
    fid.dr = min(min(0.5+vtx.x,0.5-vtx.x),min(0.5+vtx.y,0.5-vtx.y));
    if ( fid.dr < 0. ) fid.dr = 0;
    fid.dz = min(vtx.z-zextent[0],zextent[1]-vtx.z);

    break;

  default:
    MSG("EventSR",Msg::kWarning) << "Detector type " << dtype
                                 << " not supported." << endl;
    break;
  }

  return;

}

---