NtpStRecord.cxx

Go to the documentation of this file.

// 
// NtpStRecord
//
// NtpStRecord is a standard reconstruction & mc ntuple record.
//

// C++ 
#include <iostream>

// ROOT
#include "TClonesArray.h"

// MINOS headers
#include "Record/RecArrayAllocator.h"
#include "CandNtupleSR/NtpSRSlice.h"
#include "CandNtupleSR/NtpSRCluster.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRDeadChip.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "TruthHelperNtuple/NtpTHTrack.h"
#include "TruthHelperNtuple/NtpTHShower.h"
#include "MessageService/MsgService.h"

// Local package headers
#include "StandardNtuple/NtpStRecord.h"

CVSID("$Id: NtpStRecord.cxx,v 1.15 2008/02/04 02:34:29 schubert Exp $");

ClassImp(NtpStRecord)

// put using namespace std after all headers
using namespace std;

//_____________________________________________________________________________
NtpStRecord::NtpStRecord() : 
  RecRecordImp<RecCandHeader>(),vetostp(0),vetoexp(0),deadchips(0),stp(0),
   slc(0),clu(0),shw(0),trk(0),evt(0),mc(0),stdhep(0),digihit(0),thstp(0),
   thslc(0),thshw(0),thtrk(0),thevt(0) {
   // Purpose: Default constructor

  MSG("NtpSt",Msg::kVerbose) << "NtpStRecord def ctor at " << this << endl;
  this -> Init();  
}

//_____________________________________________________________________________
NtpStRecord::NtpStRecord(const RecCandHeader& hdr) : 
  RecRecordImp<RecCandHeader>(hdr),vetostp(0),vetoexp(0),deadchips(0),stp(0),slc(0),clu(0),
  shw(0),trk(0),evt(0),mc(0),stdhep(0),digihit(0),thstp(0),thslc(0),thshw(0),
  thtrk(0),thevt(0) {
  // Purpose: Normal constructor

  MSG("NtpSt",Msg::kVerbose) << "NtpStRecord normal ctor at " << this << endl;
  this -> Init(); 
}

//_____________________________________________________________________________
NtpStRecord::~NtpStRecord() {
  // Purpose: Destructor

  MSG("NtpSt",Msg::kVerbose) << "NtpStRecord dtor at " << this << endl;
  // Release arrays back to TClonesArray pool for reuse
  // Allocated memory of stored objects is retrieved via object Clear call
  RecArrayAllocator& allocator = RecArrayAllocator::Instance();
  if ( vetostp ) { allocator.ReleaseArray(vetostp); vetostp = 0; }
  if ( vetoexp ) { allocator.ReleaseArray(vetoexp); vetoexp = 0; }
  if ( deadchips ) { allocator.ReleaseArray(deadchips); deadchips = 0; }
  if ( stp ) { allocator.ReleaseArray(stp); stp = 0; }
  if ( slc ) { allocator.ReleaseArray(slc); slc = 0; }
  if ( clu ) { allocator.ReleaseArray(clu); clu = 0; }
  if ( shw ) { allocator.ReleaseArray(shw); shw = 0; }
  if ( trk ) { allocator.ReleaseArray(trk); trk = 0; }
  if ( evt ) { allocator.ReleaseArray(evt); evt = 0; }
  if ( mc ) { allocator.ReleaseArray(mc); mc = 0; }
  if ( stdhep ) { allocator.ReleaseArray(stdhep); stdhep = 0; }
  if ( digihit )  { allocator.ReleaseArray(digihit); digihit = 0; }
  if ( thstp ) { allocator.ReleaseArray(thstp); thstp = 0; }
  if ( thslc ) { allocator.ReleaseArray(thslc); thslc = 0; }
  if ( thshw ) { allocator.ReleaseArray(thshw); thshw = 0; }
  if ( thtrk ) { allocator.ReleaseArray(thtrk); thtrk = 0; }
  if ( thevt ) { allocator.ReleaseArray(thevt); thevt = 0; }

}

//_____________________________________________________________________________
void NtpStRecord::Clear(Option_t* /* option */) {
  // Purpose: Clear memory allocated to arrays so that record can
  // be reused.  

  detsim.Clear();
  if ( vetostp ) { vetostp -> Clear("C"); }
  if ( vetoexp ) { vetoexp -> Clear("C"); }
  if ( deadchips ) { deadchips->Clear("C"); }
  if ( stp ) { stp -> Clear("C"); }
  if ( slc ) { slc -> Clear("C"); }
  if ( clu ) { clu -> Clear("C"); }
  if ( shw ) { shw -> Clear("C"); }
  if ( trk ) { trk -> Clear("C"); }
  if ( evt ) { evt -> Clear("C"); }
  if ( mc ) { mc -> Clear("C"); }
  if ( stdhep ) { stdhep -> Clear("C"); }
  if ( digihit ) { digihit -> Clear("C"); }  
  if ( thstp ) { thstp -> Clear("C"); }
  if ( thslc ) { thslc -> Clear("C"); }
  if ( thshw ) { thshw -> Clear("C"); }
  if ( thtrk ) { thtrk -> Clear("C"); }
  if ( thevt ) { thevt -> Clear("C"); }

}

//_____________________________________________________________________________
void NtpStRecord::ClearStrips(Option_t* /* option */) {
  // Purpose: Clear strip and digihit information to reduce output file size

  // Standard reconstruction
  if ( vetostp ) { vetostp -> Clear("C"); } // clear array and contents
  if ( vetoexp ) { vetoexp -> Clear("C"); } //clear array and contents 
  if ( stp ) { stp -> Clear("C"); } // clear array and contents
  if ( slc ) {
    TClonesArray& slicearray = *(slc);
    for ( int islc = 0; islc <= slicearray.GetLast(); islc++ ) {
      NtpSRSlice* ntpslice = dynamic_cast<NtpSRSlice*>(slicearray.At(islc));
      ntpslice -> ClearStrips();
    }
  }

  if ( clu ) {
    TClonesArray& clusterarray = *(clu);
    for ( int iclu = 0; iclu <= clusterarray.GetLast(); iclu++ ) {
      NtpSRCluster* ntpcluster=dynamic_cast<NtpSRCluster*>
                                         (clusterarray.At(iclu));
      ntpcluster -> ClearStrips();
    }
  }

  if ( shw ) {
    TClonesArray& showerarray = *(shw);
    for ( int ishw = 0; ishw <= showerarray.GetLast(); ishw++ ) {
      NtpSRShower* ntpshower=dynamic_cast<NtpSRShower*>(showerarray.At(ishw));
      ntpshower -> ClearStrips();
    }
  }

  if ( trk ) {
    TClonesArray& trackarray = *(trk);
    for ( int itrk = 0; itrk <= trackarray.GetLast(); itrk++ ) {
      NtpSRTrack* ntptrack=dynamic_cast<NtpSRTrack*>(trackarray.At(itrk));
      ntptrack -> ClearStrips();
    }
  }

  if ( evt ) {
    TClonesArray& eventarray = *(evt);
    for ( int ievt = 0; ievt <= eventarray.GetLast(); ievt++ ) {
      NtpSREvent* ntpevent=dynamic_cast<NtpSREvent*>(eventarray.At(ievt));
      ntpevent -> ClearStrips();
    }
  }

  // Monte carlo
  if ( digihit ) { digihit -> Clear("C"); } // clear array and contents

  // Truth helper
  if ( thstp ) { thstp -> Clear("C"); } // clear array and contents

}


//_____________________________________________________________________________
void NtpStRecord::Init() {
  // 
  // Purpose: Initialize ntuple TClonesArrays
  //

  // Set variable in record base class to indicate that it is possible
  // to recover dynamic memory using Clear() method for this record type
  SetClearable(true);
  
  RecArrayAllocator& allocator = RecArrayAllocator::Instance();
  if ( !vetostp ) vetostp = allocator.GetArray("NtpSRShieldStrip");
  if ( !vetoexp ) vetoexp = allocator.GetArray("NtpSRShieldExpected"); 
  if ( !deadchips ) deadchips = allocator.GetArray("NtpSRDeadChip");
  if ( !stp ) stp = allocator.GetArray("NtpSRStrip");
  if ( !slc ) slc = allocator.GetArray("NtpSRSlice");
  if ( !clu ) clu = allocator.GetArray("NtpSRCluster");
  if ( !shw ) shw = allocator.GetArray("NtpSRShower");
  if ( !trk ) trk = allocator.GetArray("NtpSRTrack");
  if ( !evt ) evt = allocator.GetArray("NtpSREvent");
  if ( !mc )   mc = allocator.GetArray("NtpMCTruth");
  if ( !stdhep ) stdhep = allocator.GetArray("NtpMCStdHep");
  if ( !digihit ) digihit = allocator.GetArray("NtpMCDigiScintHit");
  if ( !thstp ) thstp = allocator.GetArray("NtpTHStrip");
  if ( !thslc ) thslc = allocator.GetArray("NtpTHSlice");
  if ( !thshw ) thshw = allocator.GetArray("NtpTHShower");
  if ( !thtrk ) thtrk = allocator.GetArray("NtpTHTrack");
  if ( !thevt ) thevt = allocator.GetArray("NtpTHEvent");

}

//_____________________________________________________________________________
std::ostream& NtpStRecord::Print(std::ostream& os) const {
  //
  // Purpose: Print status of ntuple record on ostream
  //

  os << "NtpStRecord::Print" << endl;
  RecRecordImp<RecCandHeader>::Print(os);

  return os;

}

//_____________________________________________________________________________
void NtpStRecord::Print(Option_t* /* option */) const {
  //
  // Purpose: Print record in form supported by TObject::Print
  //

  Print(std::cout);
  return;

}

//_____________________________________________________________________________
ReleaseType::Release_t NtpStRecord::GetRelease() const {
  //
  // Purpose: Return the Official ReleaseType determination of this record
  //

  // First try using "modern" method of retrieving release type from 
  // job history.  If that fails, try using old method to
  // retrieve release type by title.
  ReleaseType::Release_t release 
             = GetJobHistory().GetProdReleaseType(RecJobHistory::kGMinos |
                                                  RecJobHistory::kCand);
  if (release != ReleaseType::kUnknown ) return release;
    
  std::string relName = GetTitle();

  std::string mcinfo = "";
  if(GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC){
     mcinfo = "Carrot";
     string temp = mchdr.geninfo.codename;
     if(temp.size() != 0){   mcinfo = temp;  }
  }
  return ReleaseType::MakeReleaseType(relName, mcinfo);

}

//_____________________________________________________________________________
vector<const NtpSREvent *> NtpStRecord::GetEvents() const {
  // Get all events in a record. 

  vector<const NtpSREvent *> events; // empty list of events

  const int nevent = this -> evt -> GetEntriesFast();
  if (nevent < 1) return events; // return empty list if no events

   for(int i = 0; i < nevent; ++i) {
      const NtpSREvent* ntp_event 
                = dynamic_cast<const NtpSREvent*>(this -> evt -> At(i));
      if(!ntp_event) {
         MSG("NtpSt", Msg::kWarning) << "Failed to get NtpSREvent pointer at "
                                     << i << " index." << endl;
         continue;
      }
      events.push_back(ntp_event);
   }
   
   return events;
}

//_____________________________________________________________________________
vector<const NtpSRStrip *> NtpStRecord::GetStrips(const int event_index) const{
// Get strips associated with event_index.  
// If event_index < 0 (default) return all strips in record.  

   vector<const NtpSRStrip *> strips; // empty list of strips
   
   if( event_index >= this -> evt -> GetEntriesFast()) {
      MSG("NtpSt",Msg::kWarning) 
      << "GetStrips w/event index " << event_index
      << "\nis too large for evt array size " << this->evt->GetEntriesFast()
      << "." << endl;
      return strips; // empty
   }

   // if event index < 0, return all strips 
   if ( event_index < 0 ) {
      //loop over strips in record
      for(int i = 0; i < this -> stp -> GetEntriesFast(); ++i) {
         const NtpSRStrip *strip  
            = dynamic_cast<const NtpSRStrip *>(this -> stp -> At(i));
         if(!strip) {
            MSG("NtpSt", Msg::kWarning) 
            << "Failed to get NtpSRStrip strip at " << event_index << endl;
            continue;
         }       
         strips.push_back(strip);
      }
      return strips;
   }

   // event_index >= 0
   const NtpSREvent* ntp_event 
     = dynamic_cast<const NtpSREvent*>(this -> evt -> At(event_index));   
   if (!ntp_event) {
      MSG("NtpSt", Msg::kWarning) << "Failed to get NtpSREvent at "
                                  << event_index << " index." << endl;
      return strips;
   }

   const Int_t *strip_index_array = ntp_event -> stp;   
   const Int_t nstrip = ntp_event -> nstrip;

   if ( nstrip < 1 || !strip_index_array ) {
      MSG("NtpSt", Msg::kDebug) 
      << "NtpStRecord does not have valid strip index array." << endl;
      return strips;
   }
   const int stp_size = this -> stp -> GetEntriesFast();

   //loop over strips in an event
   for ( int i = 0; i < nstrip; ++i ) {
     const int index = strip_index_array[i];
     if (index >= stp_size) {
        MSG("NtpSt", Msg::kWarning) << "Invalid strip index " << index 
        << " for TClonesArray *stp of size: " << stp_size << endl;
        continue;
     }
     const NtpSRStrip *strip  
      = dynamic_cast<const NtpSRStrip *>(this -> stp -> At(index));
     if ( !strip ) {
       MSG("NtpSt", Msg::kWarning) << "Failed to get NtpSRStrip at " 
                                    << index << " index" << endl;
       continue;
     }
      
     strips.push_back(strip);
   }

   return strips;
}

//_____________________________________________________________________________
vector<const NtpSRShower *> NtpStRecord::GetShowers(const int event_index)
                                                                     const {
  // Get showers associated with event_index. 
  // If event_index < 0 (default) return all showers in record.

   vector<const NtpSRShower *> showers; // empty list of showers

   if ( event_index >= this -> evt -> GetEntriesFast() ) {
      MSG("NtpSt",Msg::kWarning) 
      << "GetShowers w/event index " << event_index
      << "\nis too large for evt array size " << this->evt->GetEntriesFast()
      << "." << endl;
      return showers; // empty shower list      
   }

   if ( event_index < 0 ) {
      //loop over showers in a record
      for ( int i = 0; i < this -> shw -> GetEntriesFast(); ++i ) {
        const NtpSRShower *shower  
               = dynamic_cast<const NtpSRShower *>(this -> shw -> At(i));
        if(!shower) {
           MSG("NtpSt", Msg::kWarning) << "No NtpSRShower pointer at " 
                                       << event_index << " index." << endl;
           continue;
        }
        showers.push_back(shower);
     }
     return showers; // all showers in record
   }

   // event_index >= 0
   const NtpSREvent* ntp_event 
          = dynamic_cast<const NtpSREvent*>(this -> evt -> At(event_index));   
   if(!ntp_event) {
      MSG("NtpSt", Msg::kWarning) << "Failed to get NtpSREvent pointer at "
                                  << event_index << " index." << endl;
      return showers;
   }

   const Int_t *shower_index_array = ntp_event -> shw;   
   const Int_t nshower = ntp_event -> nshower;

   if(nshower < 1 || !shower_index_array) {
      MSG("NtpSt", Msg::kDebug) 
      << "NtpStRecord does not have valid shower index array." <<endl;
      return showers;
   }

   const int shw_size = this -> shw -> GetEntriesFast();

   //loop over showers in an event
   for(int i = 0; i < nshower; ++i) {
      const int index = shower_index_array[i];
      if(index >= shw_size) {
         MSG("NtpSt", Msg::kWarning) << "Invalid shower index " << index 
         << " for TClonesArray *shw of size: " << shw_size << endl;
         continue;
      }
      const NtpSRShower *shower  
               = dynamic_cast<const NtpSRShower *>(this -> shw -> At(index));
      if(!shower) {
         MSG("NtpSt", Msg::kWarning) << "No NtpSRShower pointer at " 
                                     << index << " index." << endl;
         continue;
      }
      
      showers.push_back(shower);
   }

   return showers;
}

//_____________________________________________________________________________
vector<const NtpSRTrack *> NtpStRecord::GetTracks(const int event_index) const{
  // Get tracks associated with event_index. 
  // If event_index < 0 return all tracks in record.

   vector<const NtpSRTrack *> tracks; // empty list of tracks

   if(event_index >= this -> evt -> GetEntriesFast()) {
      MSG("NtpSt",Msg::kWarning) 
      << "GetTracks w/event index " << event_index
      << "\nis too large for evt array size " << this -> evt->GetEntriesFast()
      << "." << endl;
      return tracks;
   }

   if( event_index < 0 ) {
     //loop over tracks in a record
     for(int i = 0; i < this -> trk -> GetEntriesFast(); ++i) {
       const NtpSRTrack *track  
                = dynamic_cast<const NtpSRTrack *>(this -> trk -> At(i));
       if (!track) {
            MSG("NtpSt", Msg::kWarning) << "No track at " << event_index 
            << endl;
            continue;
       }
         
       tracks.push_back(track);
     }
     return tracks; // all tracks in record
   }
   
   // event_index >= 0
   const NtpSREvent* ntp_event 
            = dynamic_cast<const NtpSREvent*>(this -> evt -> At(event_index));
   if(!ntp_event) {
     MSG("NtpSt", Msg::kWarning) << "Failed to find NtpSREvent pointer at "
                                  << event_index << " index." << endl;
      return tracks;
   }

   const Int_t *track_index_array = ntp_event -> trk;
   const Int_t ntrack = ntp_event -> ntrack;

   if(ntrack < 1 || !track_index_array) {
      MSG("NtpSt", Msg::kDebug) 
      << "NtpStRecord does not have valid track index array." <<endl;
      return tracks;
   }

   const int trk_size = this -> trk -> GetEntriesFast();

   //loop over tracks in an event
   for(int i = 0; i < ntrack; ++i) {
      const int index = track_index_array[i];
      if ( index >= trk_size ) {
         MSG("NtpSt", Msg::kWarning) << "Invalid track index " << index 
         << " for TClonesArray *trk of size: " << trk_size << endl;
         continue;
      }
      const NtpSRTrack *track  
          = dynamic_cast<const NtpSRTrack *>(this -> trk -> At(index));
      if(!track) {
         MSG("NtpSt", Msg::kWarning) 
         << "Failed to get NtpSRTrack pointer at " << index << endl;
         continue;
      }
      
      tracks.push_back(track);
   }

   return tracks;
}

//_____________________________________________________________________________
const NtpMCTruth* NtpStRecord::GetEventMCTruth(const int event_index) const {
  // Get NtpMCTruth entry corresponding to truth helper identified best
  // match for specified event_index.

   if(event_index < 0 || event_index >= this -> thevt -> GetEntriesFast()) {
      // 1:1 relationship between thevt and evt entries, so treat this
      // error seriously  
      MSG("NtpSt",Msg::kWarning) 
      << "GetEventMCTruth w/event index " << event_index
      << "\nis too large for thevt array size " 
      << this -> thevt->GetEntriesFast() << "." << endl;
      return 0;
   }
   
   const NtpTHEvent* ntp_event 
      = dynamic_cast<const NtpTHEvent*>(this -> thevt -> At(event_index));
   if(!ntp_event) {
      MSG("NtpSt", Msg::kWarning) << "Failed to get NtpTHEvent pointer at "
                                  << event_index << " index." << endl;
      return 0;
   }

   const int mctruth_index = ntp_event -> neumc;
   if ( mctruth_index < 0 ) return 0;  // valid response if no best match

   if(  mctruth_index >= this -> mc -> GetEntriesFast() ) {
      MSG("NtpSt", Msg::kError) 
      << "Invalid thevt.neumc index " << mctruth_index 
      << " exceeds bounds of " 
      << "mc array, size " << this -> mc -> GetEntriesFast() << endl;
      return 0;
   }
   
   const NtpMCTruth *truth  
          = dynamic_cast<const NtpMCTruth *>(this -> mc -> At(mctruth_index));
   if(!truth) {
      MSG("NtpSt", Msg::kDebug) << "No mc truth at " << mctruth_index << endl;
   }
   
   return truth;
}

//_____________________________________________________________________________
vector<const NtpMCTruth *> NtpStRecord::GetMCTruths() const {
  // Get all mc truth events in the record. 

  vector< const NtpMCTruth* > mctruths; // empty list of mctruths

  // Push ptrs to all mc truth entries to output vector 
  for ( int imc = 0; imc < mc->GetEntriesFast(); imc++ ) {
    const NtpMCTruth* ntpmctruth 
                = dynamic_cast<const NtpMCTruth*>( mc -> At(imc) );
    mctruths.push_back(ntpmctruth);
  }
   
  return mctruths;

}

//_____________________________________________________________________________
vector<const NtpMCStdHep *> NtpStRecord::GetMCStdHeps(int mctruth_index) const{
  // Get stdhep entries associated with mc truth entry mctruth_index. 
  // If mctruth_index < 0 (default) return all stdhep entries in record.

  vector<const NtpMCStdHep *> stdheps; // empty list of stdhep entries

  // Determine stdhep indices of interest
  Int_t begstdhep_index = -1;
  Int_t endstdhep_index = -1;
  if ( mctruth_index < 0 ) {
    // If mc truth entry has not been specified, sthep indices of interest
    // is entire range
    begstdhep_index = 0;
    endstdhep_index = (stdhep -> GetEntriesFast()) - 1;
  }
  else {
    // stdhep indices of interest are those associated with mctruth_index
    if ( mctruth_index >= mc -> GetEntriesFast() ) {
      MSG("NtpSt",Msg::kWarning) 
             << "GetMCStdHeps w/mc truth array entry " << mctruth_index
             << "\nis too large for mc array size " 
             << mc -> GetEntriesFast()
             << "." << endl;
      return stdheps; // empty list of stdhep entries
    }
    const NtpMCTruth* ntpmctruth 
          = dynamic_cast<const NtpMCTruth*>(mc -> At(mctruth_index));
    if ( !ntpmctruth ) {
      MSG("NtpSt", Msg::kWarning) << "Failed to find NtpMCTruth pointer at "
                                  << mctruth_index << " index." << endl;
      return stdheps; // empty list of stdhep entries
    }

    // set stdhep indices according to this mc entry
    begstdhep_index = ntpmctruth->stdhep[0];
    endstdhep_index = ntpmctruth->stdhep[1];
  }
  
  if ( begstdhep_index < 0 || endstdhep_index < 0 ) {
    // empty stdheps list, don't think this should ever happen legitimately
    MSG("NtpSt",Msg::kWarning) << "No stdhep entries associated with "
                               << "mctruth_index " << mctruth_index << endl;
    return stdheps;
  }  
 
  // Loop over selected stdhep indices and push to output vector 
  for ( int istd = begstdhep_index; istd <= endstdhep_index; istd++ ) {
    const NtpMCStdHep *ntpstdhep  
          = dynamic_cast<const NtpMCStdHep *>( stdhep -> At(istd) );
    if ( !ntpstdhep ) {
      MSG("NtpSt", Msg::kWarning) << "No stdhep at index " << istd << endl;
      continue;
    }
    stdheps.push_back(ntpstdhep);
  }

  return stdheps; // stdhep entries in selected index range

}

//_____________________________________________________________________________
vector< const NtpTHEvent* > NtpStRecord::GetTHEvents() const {
  // Get all mc truth helper events in the record. 

  vector< const NtpTHEvent* > thevents; 

  // Fill ptrs to all truth helper events to the output vector
  for ( int ith = 0; ith < thevt->GetEntriesFast(); ith++ ) {
    const NtpTHEvent* ntpthevent 
                = dynamic_cast<const NtpTHEvent*>( thevt -> At(ith) );
    thevents.push_back(ntpthevent);
  }
   
  return thevents;

}

//_____________________________________________________________________________
vector< const NtpTHTrack* > NtpStRecord::GetTHTracks( int event_index ) const{
  // Get truth helper tracks associated with event_index. 
  // If event_index < 0 return all truth helper tracks in record.

  vector<const NtpTHTrack *> thtracks; // empty list of tracks

  if ( thtrk->GetEntriesFast() <= 0 ) return thtracks; // no th tracks

  if ( event_index < 0 ) {
    // Fill all th tracks to output vector
    for ( int ith = 0; ith < thtrk -> GetEntriesFast(); ith++ ) {
      const NtpTHTrack *thtrack  
            = dynamic_cast<const NtpTHTrack *>(thtrk -> At(ith));
      thtracks.push_back(thtrack);
    }
  }
  else {
    if ( event_index >= evt -> GetEntriesFast() ) {
      MSG("NtpSt",Msg::kWarning) 
      << "GetTHTracks w/event index " << event_index
      << "\nis too large for evt array size " << evt -> GetEntriesFast()
      << "." << endl;
      return thtracks; // empty list
    }
 
    std::vector<const NtpSRTrack*> tracks = GetTracks(event_index);
    for ( unsigned int itrk = 0; itrk < tracks.size(); itrk++ ) {
      const NtpSRTrack* ntptrack = tracks[itrk];
      Int_t trkindex = ntptrack->index;
      const NtpTHTrack* ntpthtrack = dynamic_cast<const NtpTHTrack*>
                                     (thtrk->At(trkindex));
      thtracks.push_back(ntpthtrack);
    }
  }

  return thtracks; // return filled tracks

}

//_____________________________________________________________________________
vector< const NtpTHShower* > NtpStRecord::GetTHShowers( int event_index) const{
  // Get truth helper showers associated with event_index. 
  // If event_index < 0 return all truth helper showers in record.

  vector<const NtpTHShower* > thshowers; // empty list of showers

  if ( thshw->GetEntriesFast() <= 0 ) return thshowers; // no th showers

  if ( event_index < 0 ) {
    // Fill all th showers to output vector
    for ( int ishw = 0; ishw < thshw -> GetEntriesFast(); ishw++ ) {
      const NtpTHShower *thshower  
            = dynamic_cast<const NtpTHShower* >(thshw -> At(ishw));
      thshowers.push_back(thshower);
    }
  }
  else {
    if ( event_index >= evt -> GetEntriesFast() ) {
      MSG("NtpSt",Msg::kWarning) 
      << "GetTHShowers w/event index " << event_index
      << "\nis too large for evt array size " << evt -> GetEntriesFast()
      << "." << endl;
      return thshowers; // empty list
    }
 
    std::vector<const NtpSRShower*> showers = GetShowers(event_index);
    for ( unsigned int ishw = 0; ishw < showers.size(); ishw++ ) {
      const NtpSRShower* ntpshower = showers[ishw];
      Int_t shwindex = ntpshower->index;
      const NtpTHShower* ntpthshower = dynamic_cast<const NtpTHShower*>
                                     (thshw->At(shwindex));
      thshowers.push_back(ntpthshower);
    }
  }

  return thshowers; // return filled showers

}

---