MadBase.cxx

Go to the documentation of this file.

#ifndef madbase_cxx
#define madbase_cxx
#include <iostream>
#include <cmath>
//root includes
#include "TClonesArray.h"
#include "Mad/MadBase.h"
#include "Conventions/Munits.h"
#include "DataUtil/EnergyCorrections.h"
//ClassImp(MadBase)

using namespace EnergyCorrections;

MadBase::MadBase(TChain *chainSR,TChain *chainMC,
                 TChain *chainTH,TChain *chainEM)
{

  fCurRun = -1;
  lastEntry = -1;
  if(!chainSR) {
    record = 0;
    emrecord = 0;
    mcrecord = 0;
    threcord = 0;
    strecord = 0;
    Clear();
    whichSource = -1;
    isMC = true;
    isTH = true;
    isEM = true;
    isST = false;
    Nentries = -1;
    fCurRun = -1;
    return;
  }

  InitChain(chainSR,chainMC,chainTH,chainEM);
  whichSource = 0;

}

MadBase::MadBase(JobC *j,string path,int entries)
{
  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  strecord = 0;
  Clear();
  isMC = true;
  isTH = true;
  isEM = true;
  isST = false;
  Nentries = entries;
  whichSource = 1;
  jcPath = path;
  fCurRun = -1;
  lastEntry = -1;
  fJC = j;
  fChain = NULL;
}

MadBase::~MadBase()
{
  if(fChain) delete fChain;
}

Int_t MadBase::GetEntry(Int_t run,Int_t snarl)
{
  if(whichSource==0) return GetEntry(fChain->GetEntryNumber(run,snarl));
  return 0;
}

Int_t MadBase::GetEntry(Int_t entry)
{
  //zero everything before next snarl
  Clear();
  
  //get records
  int status = 0;
  if(whichSource==-1) return 0;
  else if(whichSource==0) {
    if(entry!=lastEntry) {
      status = fChain->Get(entry);
      lastEntry = entry;
    }
    else status = 1;
    record = fChain->Record;
    strecord = fChain->stRecord;
    mcrecord = fChain->mcRecord;
    threcord = fChain->thRecord;
    emrecord = fChain->emRecord;
  }
  else if(whichSource==1) {
    //allow for possibility that reco path can change event by event:
    isMC = isTH = isEM = true; 
    //Get entry:
    //cout << "Getting Cand entry " << entry << " from " 
    // << jcPath.c_str() << endl;
    if(entry!=lastEntry) {
      if(entry<lastEntry) fJC->Input.Prev(lastEntry-entry-2);
      //else fJC->Input.Next(entry-lastEntry);
      fJC->Path(jcPath.c_str()).RunNin(1);
      lastEntry = entry;
    }
    //Look for records in mom:
    record   = (NtpSRRecord*) fJC->Mom.GetFragment("NtpSRRecord");
    mcrecord = (NtpMCRecord*) fJC->Mom.GetFragment("NtpMCRecord");
    threcord = (NtpTHRecord*) fJC->Mom.GetFragment("NtpTHRecord");
    emrecord = (NtpEMRecord*) fJC->Mom.GetFragment("NtpEMRecord");
    //if no SR record quit
    if(record==0) return 0;
    //otherwise check for other records:
    status = 1;
    if(mcrecord==0) isMC = false;
    if(threcord==0) isTH = false;
    if(emrecord==0) isEM = false;
  }
  else return 0;
  
  //set other things:
  if(isST){
    ntpHeader = &(strecord->GetHeader());
    mcHeader = &(strecord->mchdr);
    shieldSummary = &(strecord->vetohdr);
    eventSummary = &(strecord->evthdr);
    dmxStatus = &(strecord->dmxstatus);
    detStatus = &(strecord->detstatus);
    if(mcHeader->nmc==0) {
      isMC = false; isTH = false;
    }
    else {
      isMC = true;
      TClonesArray &array = *(strecord->thstp);
      if(array.GetEntries()>0) isTH = true;      
    }
  }
  else {
    ntpHeader = &(record->GetHeader());
    shieldSummary = &(record->vetohdr);
    eventSummary = &(record->evthdr);
    dmxStatus = &(record->dmxstatus);
    detStatus = &(record->detstatus);
  }
  fCurRun = ntpHeader->GetRun();

  LoadSlice(0);  
  LoadEvent(0);
  LoadTrack(0);
  LoadCluster(0);
  LoadShower(0);
  LoadShieldStrip(0);

  if(isEM) {
    ntpEMSummary = &(emrecord->emhdr);
    LoadEMShower(0);
  }

  if(isMC) {
    if(isST) mcHeader = &(strecord->mchdr);
    else mcHeader = &(mcrecord->mchdr);
    LoadTruth(0);
    LoadStdHep(0);
  }

  if(isTH) {
    LoadTHSlice(0);
    LoadTHEvent(0);
    LoadTHTrack(0);
    LoadTHShower(0);
  }
  return status;

}

void MadBase::InitChain(TChain *chainSR,TChain *chainMC,
                        TChain *chainTH,TChain *chainEM, bool build_index)
{

  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  strecord = 0;

  fChain = new MadChain(chainSR,chainMC,chainTH,chainEM,build_index);
  isMC = fChain->isMC;
  isTH = fChain->isTH;
  isEM = fChain->isEM;
  isST = fChain->isST;

  Nentries = fChain->Nentries;  

  if(isST) std::cout << "ST Present" << std::endl;
  if(isMC) std::cout << "MC Present" << std::endl;
  if(isTH) std::cout << "TH Present" << std::endl;
  if(isEM) std::cout << "EM Present" << std::endl;
  
}

void MadBase::Clear(){

  ntpHeader = 0;
  shieldSummary = 0;
  eventSummary = 0;
  ntpSlice=0;
  ntpEvent=0;
  ntpTrack=0;
  ntpCluster=0;
  ntpShower=0;
  ntpStrip=0;
  ntpShieldStrip=0;

  ntpTruth=0;
  ntpStdHep=0;

  ntpTHSlice=0;
  ntpTHTrack=0;
  ntpTHShower=0;

  ntpEMSummary = 0;
  ntpEMShower = 0;

}


Bool_t MadBase::LoadTrack(Int_t itrk){

  Bool_t status = false;
  TClonesArray* pointTrackArray = NULL;
  if(isST) pointTrackArray = (strecord->trk);
  else pointTrackArray = (record->trk);
  TClonesArray& trackArray = *pointTrackArray;
//MAK:06/06/05  if(itrk>=0&&itrk<trackArray.GetEntries()) {
  if(itrk>=0&&itrk<trackArray.GetEntriesFast()) {
    ntpTrack = dynamic_cast<NtpSRTrack *>(trackArray[itrk]);
    status = true;
  }
  else ntpTrack = 0;
  
  return status;
}

Bool_t MadBase::LoadLargestTrackFromEvent(Int_t ievt,Int_t &itrk){

  Bool_t status = false;  
  ntpTrack = 0;
  if(!LoadEvent(ievt)) return status;

  Float_t largest_ph = 0;
  int *tracks = ntpEvent->trk;  
  for(int i=0;i<ntpEvent->ntrack;i++){    
    if(!LoadTrack(tracks[i])) continue;
    // MAK: trklen should be an absolute value
    // MAK: Also, it shouldn't matter for beam events
    Int_t trklen=abs(ntpTrack->plane.end-ntpTrack->plane.beg+1);
    if(trklen>largest_ph) {
      itrk = tracks[i];
      largest_ph = trklen;
    }
  }
  ntpTrack = 0;
  if(largest_ph>0){
    status = true;
    LoadTrack(itrk);
  } 
  return status;
}

Bool_t MadBase::LoadTHEvent(Int_t ievt){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray* pointEventArray = NULL;
  if(isST) pointEventArray = (strecord->thevt);
  else pointEventArray = (threcord->thevt);
  TClonesArray& eventArray = *pointEventArray;
  //MAK:06/06/05  if(ievt>=0&&ievt<eventArray.GetEntries()) {
  if(ievt>=0&&ievt<eventArray.GetEntriesFast()) {
    ntpTHEvent = dynamic_cast<NtpTHEvent *>(eventArray[ievt]);
    status = true;
  }
  else ntpTHEvent = 0;
  
  return status;
}

Bool_t MadBase::LoadTHTrack(Int_t itrk){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray* pointTrackArray = NULL;
  if(isST) pointTrackArray = (strecord->thtrk);
  else pointTrackArray = (threcord->thtrk);
  TClonesArray& trackArray = *pointTrackArray;
//MAK:06/06/05  if(itrk>=0&&itrk<trackArray.GetEntries()) {
  if(itrk>=0&&itrk<trackArray.GetEntriesFast()) {
    ntpTHTrack = dynamic_cast<NtpTHTrack *>(trackArray[itrk]);
    status = true;
  }
  else ntpTHTrack = 0;
  
  return status;
}

Bool_t MadBase::LoadCluster(Int_t iclu){
 
  Bool_t status = false;
  TClonesArray* pointClusterArray = NULL;
  if(isST) pointClusterArray = (strecord->clu);
  else pointClusterArray = (record->clu);
  TClonesArray& clusterArray = *pointClusterArray;
//MAK:06/06/05  if(iclu>=0&&iclu<clusterArray.GetEntries()) {
  if(iclu>=0&&iclu<clusterArray.GetEntriesFast()) {
    ntpCluster = dynamic_cast<NtpSRCluster *>(clusterArray[iclu]);
    status = true;
  }
  else ntpCluster = 0;
   
  return status;
}

Bool_t MadBase::LoadShower(Int_t ishw){

  Bool_t status = false;
  TClonesArray* pointShowerArray = NULL;
  if(isST) pointShowerArray = (strecord->shw);
  else pointShowerArray = (record->shw);
  TClonesArray& showerArray = *pointShowerArray;
//MAK:06/06/05  if(ishw>=0&&ishw<showerArray.GetEntries()) {
  if(ishw>=0&&ishw<showerArray.GetEntriesFast()) {
    ntpShower = dynamic_cast<NtpSRShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpShower = 0;
  
  return status;
}

Bool_t MadBase::LoadLargestShowerFromEvent(Int_t ievt,Int_t &ishw){

  Bool_t status = false;
  ntpShower = 0;
  if(!LoadEvent(ievt)) return status;
  
  Float_t largest_ph = 0;
  int *showers = ntpEvent->shw;  
  for(int i=0;i<ntpEvent->nshower;i++){    
    if(!LoadShower(showers[i])) continue;
    if(ntpShower->ph.gev>largest_ph) {
      ishw = showers[i];
      largest_ph = ntpShower->ph.gev;
    }
  }
  ntpShower = 0;
  if(largest_ph>0){
    status = true;
    LoadShower(ishw);
  } 
  return status;
}

Bool_t MadBase::LoadShower_Jim(Int_t ievt,Int_t &ishw, Int_t detector){

  Bool_t status = false;
  ntpShower = 0;
  if(!LoadEvent(ievt)) return status;
  int *showers = ntpEvent->shw;  
  if(ntpEvent->ntrack==0 && ntpEvent->nshower>0) {
   ishw = showers[0];
   ntpShower = 0;  
   status = true;
   LoadShower(ishw);
  }
  if(ntpEvent->primshw<0) return false;

  if(ntpEvent->nshower>0 && ntpEvent->ntrack> 0) {
   if(!LoadShower(showers[0])) return status;
   ishw = showers[0];
   ntpShower = 0;    
   LoadShower(ishw);
   int track_index   = -1;
   LoadLargestTrackFromEvent(ievt,track_index);  
   Double_t reco_eshwcc,reco_eshwccw; 
   Int_t shwtype=0;  
   shwtype=0;
   reco_eshwcc       = RecoShwEnergy(ishw,shwtype,detector);
   shwtype=1;
   reco_eshwccw      = RecoShwEnergy(ishw,shwtype,detector); 
   if(track_index!=-1){
      Double_t trvz= ntpTrack ->vtx.z; 
      Double_t shvz= ntpShower->vtx.z;
      if(fabs(trvz-shvz)<0.5 || (fabs(trvz-shvz)>=0.5 && reco_eshwcc > 2)) {
       status=true;       
      }
      else if(fabs(trvz-shvz)>=0.5 && reco_eshwcc <= 2) {
       status=false;      
     }
    }
  }
  if(!status) ishw=-1;
  
  return status;
}
//
Float_t MadBase::RecoShwEnergy(Int_t ishw, Int_t opt, Int_t det){
  //use SR reco
  Float_t shower_ene=0;
                                                                                                                                         
  if(LoadShower(ishw)) {
  if(det==1){
    if(opt==-1) shower_ene = ntpShower->ph.gev;
    if(opt==0)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linCCgev,CandShowerHandle::kCC);
    if(opt==1)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC);
    if(opt==2)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linNCgev,CandShowerHandle::kNC);
    if(opt==3)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC);
    if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
    return shower_ene;
  }
   if(det==2){
    if(opt==-1) shower_ene = ntpShower->ph.gev;
    if(opt==0)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linCCgev,CandShowerHandle::kCC);
    if(opt==1)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC);
    if(opt==2)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linNCgev,CandShowerHandle::kNC);
    if(opt==3)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC);
    if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
    return shower_ene;
                                                                                                                                         
  }
}
  return 0;
}


//
Bool_t MadBase::LoadShowerAtTrackVertex(Int_t ievt,Int_t itrk,Int_t &ishw){

  Bool_t status = false;
  ntpShower = 0;
  if(!LoadEvent(ievt)) return status;
  if(!LoadTrack(itrk)) return status;

  Float_t closest_ph = 0;
  Float_t closest_vtx = 1000.;
  const Float_t close_enough=7*0.06; // about 1 interaction length
  int *showers = ntpEvent->shw;
  for(int i=0;i<ntpEvent->nshower;i++){
    if(!LoadShower(showers[i])) continue;
    /*
    Float_t vtx_sep = TMath::Sqrt(TMath::Power(ntpShower->vtx.x-ntpTrack->vtx.x,2) +
                                  TMath::Power(ntpShower->vtx.y-ntpTrack->vtx.y,2) +
                                  TMath::Power(ntpShower->vtx.z-ntpTrack->vtx.z,2));
    */
    Float_t vtx_sep = TMath::Abs(ntpShower->vtx.z-ntpTrack->vtx.z);

    // if this is the closest shower
    if(vtx_sep < closest_vtx) {
      ishw = showers[i];
      closest_vtx = vtx_sep;
      // this is a check to make sure we have R1.16 or greater shw variables
      if(ntpShower->shwph.linCCgev>0) closest_ph = ntpShower->shwph.linCCgev;
      else closest_ph=ntpShower->ph.gev;
    }
    // if this isn't the closest shower but it's within close_enough m of the
    // track vertex and has a larger pulseheight
    else if(vtx_sep<close_enough && ntpShower->ph.gev>closest_ph){
      ishw = showers[i];
      closest_vtx = vtx_sep;

      if(ntpShower->shwph.linCCgev>0) closest_ph = ntpShower->shwph.linCCgev;
      else closest_ph=ntpShower->ph.gev;
    }
  }
  ntpShower = 0;
  // only associate showers within close_enough m of the track vertex
  if(closest_vtx<close_enough){
    status = true;
    LoadShower(ishw);
  }
  return status;
}

Bool_t MadBase::LoadTHShower(Int_t ishw){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray* pointShowerArray = NULL;
  if(isST) pointShowerArray = (strecord->thshw);
  else pointShowerArray = (threcord->thshw);
  TClonesArray& showerArray = *pointShowerArray;
//MAK:06/06/05  if(ishw>=0&&ishw<showerArray.GetEntries()) {
  if(ishw>=0&&ishw<showerArray.GetEntriesFast()) {
    ntpTHShower = dynamic_cast<NtpTHShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpTHShower = 0;
  
  return status;
}

Bool_t MadBase::LoadSlice(Int_t islc){

  Bool_t status = false;
  TClonesArray* pointSliceArray = NULL;
  if(isST) pointSliceArray = (strecord->slc);
  else pointSliceArray = (record->slc);
  TClonesArray& sliceArray = *pointSliceArray;
//MAK:06/06/05  if(islc>=0&&islc<sliceArray.GetEntries()) {
  if(islc>=0&&islc<sliceArray.GetEntriesFast()) {
    ntpSlice = dynamic_cast<NtpSRSlice *>(sliceArray[islc]);
    status = true;
  }
  else ntpSlice = 0;
  
  return status;
}

Bool_t MadBase::LoadTHSlice(Int_t islc){

  Bool_t status = false;
  if(!isTH) return status;
  TClonesArray* pointSliceArray = NULL;
  if(isST) pointSliceArray = (strecord->thslc);
  else pointSliceArray = (threcord->thslc);
  TClonesArray& sliceArray = *pointSliceArray;
//MAK:06/06/05  if(islc>=0&&islc<sliceArray.GetEntries()) {
  if(islc>=0&&islc<sliceArray.GetEntriesFast()) {
    ntpTHSlice = dynamic_cast<NtpTHSlice *>(sliceArray[islc]);
    status = true;
  }
  else ntpTHSlice = 0;
  
  return status;
}

Bool_t MadBase::LoadEvent(Int_t ievt){

  Bool_t status = false;
  TClonesArray* pointEventArray = NULL;
  if(isST) pointEventArray = (strecord->evt);
  else pointEventArray = (record->evt);
  TClonesArray& eventArray = *pointEventArray;
//MAK:06/06/05  if(ievt>=0&&ievt<eventArray.GetEntries()) {
  if(ievt>=0&&ievt<eventArray.GetEntriesFast()) {
    ntpEvent = dynamic_cast<NtpSREvent *>(eventArray[ievt]);
    status = true;
  }
  else ntpEvent = 0;
  
  return status;
}

Bool_t MadBase::LoadTruth(Int_t itr){
  
  Bool_t status = false;
  if(!isMC) return status;
  TClonesArray* pointMcArray = NULL;
  if(isST) pointMcArray = (strecord->mc);
  else pointMcArray = (mcrecord->mc);
  TClonesArray& mcArray = *pointMcArray;
//MAK:06/06/05  if(itr>=0&&itr<mcArray.GetEntries()) {
  if(itr>=0&&itr<mcArray.GetEntriesFast()) {
    ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[itr]);
    status = true;
  }
  else ntpTruth = 0;
  
  return status;
}

Bool_t MadBase::LoadStdHep(Int_t ish){
  
  Bool_t status = false;
  if(!isMC) return status;
  TClonesArray* pointStdhepArray = NULL;
  if(isST) pointStdhepArray = (strecord->stdhep);
  else pointStdhepArray = (mcrecord->stdhep);
  TClonesArray& stdhepArray = *pointStdhepArray;
//MAK:06/06/05  if(ish>=0&&ish<stdhepArray.GetEntries()) {
  if(ish>=0&&ish<stdhepArray.GetEntriesFast()) {
    ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[ish]);
    status = true;
  }
  else ntpStdHep = 0;
  
  return status;
}

Bool_t MadBase::LoadTruthForReco(Int_t ievt,Int_t &itruth){

  Bool_t status = false;  
  ntpTruth = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;

  TClonesArray* pointMcArray = NULL;
  if(isST) pointMcArray = (strecord->mc);
  else pointMcArray = (mcrecord->mc);
  TClonesArray& mcArray = *pointMcArray;
  if(mcArray.GetEntries()==1) { //if only one MC event is present just load it
    LoadTruth(0);
    return true;
  }

  float mindist = 999;
  int itr = 0;
  //loop over all MC events and look for closest vertex to reco vertex
  for(int i=0;i<mcArray.GetEntries();i++){
    ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[i]);
    float dist = sqrt(TMath::Power(ntpTruth->vtxx/100.-ntpEvent->vtx.x,2)
                      +TMath::Power(ntpTruth->vtxy/100.-ntpEvent->vtx.y,2)
                      +TMath::Power(ntpTruth->vtxz/100.-ntpEvent->vtx.z,2));
    if(dist<mindist) { 
      itr = i;
      mindist = dist;
    }
  }

  if(mindist<0.5) {  //Need to place a fairly loose cut to pick up showers
    status = true;   //From looking at FarDet events, a cut of 0.5m 
    LoadTruth(itr);  //correctly pairs up 95% of events (LE beam,unoscillated) 
    itruth = itr;
  }
  return status;
}

Bool_t MadBase::LoadSliceForRecoTH(Int_t ievt,Int_t &islice){

  Bool_t status = false;
  ntpSlice = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;
  if(!isTH) return status;

  /*
    Int_t itru = 0;
    if(LoadTruthForRecoTH(ievt,itru)){
    Int_t count = 0;
    while(!status){
    if(!LoadTHSlice(count)) break;
    if(ntpTHSlice->neumc==itru) {
    if(LoadSlice(count)){
    islice = count;
    status = true;
    }
    }
    count += 1;
    }
    }
  */

  //Tingjun's sugggested change:
  islice = ntpEvent->slc;
  if (LoadSlice(islice)&&LoadTHSlice(islice)) status = true;
  
  return status;
}


Bool_t MadBase::LoadTruthForRecoTH(Int_t ievt,Int_t &itruth){

  Bool_t status = false;
  ntpTruth = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;
  if(!isTH) return status;
  
  if(ntpEvent->ntrack>0) {
    int *tracks = ntpEvent->trk;
    if(LoadTHTrack(int(tracks[0]))){
      int ind = ntpTHTrack->neumc;
      if(LoadTruth(ind)) {
        itruth = ind;
        status = true;
      }
    }
  }
  if(!status&&ntpEvent->nshower>0) {
    int *showers = ntpEvent->shw;    
    if(LoadTHShower(int(showers[0]))){
      int ind = ntpTHShower->neumc;
      if(LoadTruth(ind)) {
        itruth = ind;
        status = true;
      }
    }
  }
  return status;

}

Bool_t MadBase::LoadStrip(Int_t istp){

  Bool_t status = false;
  TClonesArray* pointStripArray = NULL;
  if(isST) pointStripArray = (strecord->stp);
  else pointStripArray = (record->stp);
  TClonesArray& stripArray = *pointStripArray;
  //MAK:06/06/05 if(istp>=0&&istp<stripArray.GetEntries()) { 
  if(istp>=0&&istp<stripArray.GetEntriesFast()) {
    ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[istp]);
    status = true;
  }
  else ntpStrip = 0;
  
  return status;
}

Bool_t MadBase::LoadShieldStrip(Int_t istp){

  Bool_t status = false;
  TClonesArray* pointStripArray = NULL;
  if(isST) pointStripArray = (strecord->vetostp);
  else pointStripArray = (record->vetostp);
  TClonesArray& stripArray = *pointStripArray;
//MAK:06/06/05  if(istp>=0&&istp<stripArray.GetEntries()) {
  if(istp>=0&&istp<stripArray.GetEntriesFast()) {
    ntpShieldStrip = dynamic_cast<NtpSRShieldStrip *>(stripArray[istp]);
    status = true;
  }
  else ntpShieldStrip = 0;
  
  return status;
}

Bool_t MadBase::LoadEMShower(Int_t ishw){

  Bool_t status = false;
  if(!isEM) return status; 
  TClonesArray& showerArray = *(emrecord->emshw);
//MAK:06/06/05  if(ishw>=0&&ishw<showerArray.GetEntries()) {
  if(ishw>=0&&ishw<showerArray.GetEntries()) {
    ntpEMShower = dynamic_cast<NtpEMShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpEMShower = 0;
  
  return status;
}

const NtpSRTrack* MadBase::GetTrack(Int_t i){
  if(!LoadTrack(i)) return 0;
  else return ntpTrack;

}

const NtpSRTrack* MadBase::GetLargestTrackFromEvent(Int_t i,Int_t& j){
  if(!LoadLargestTrackFromEvent(i,j)) return 0;
  else return ntpTrack;
}
const NtpSRShower* MadBase::GetShower(Int_t i){
  if(!LoadShower(i)) return 0;
  else return ntpShower;
}
const NtpSRCluster* MadBase::GetCluster(Int_t i){
  if(!LoadCluster(i)) return 0;
  else return ntpCluster;
}
const NtpSRShower* MadBase::GetLargestShowerFromEvent(Int_t i,Int_t& j){
  if(!LoadLargestShowerFromEvent(i,j)) return 0;
  else return ntpShower;
}
const NtpSRShower* MadBase::GetShowerAtTrackVertex(Int_t i,Int_t j,Int_t& k){
  if(!LoadShowerAtTrackVertex(i,j,k)) return 0;
  else return ntpShower;
}
const NtpSREvent* MadBase::GetEvent(Int_t i){
  if(!LoadEvent(i)) return 0;
  else return ntpEvent;
}
const NtpSRSlice* MadBase::GetSlice(Int_t i){
  if(!LoadSlice(i)) return 0;
  else return ntpSlice;
}
const NtpSRShieldStrip* MadBase::GetShieldStrip(Int_t i){
  if(!LoadShieldStrip(i)) return 0;
  else return ntpShieldStrip;
  
}
const NtpSRStrip* MadBase::GetStrip(Int_t i){
  if(!LoadStrip(i)) return 0;
  else return ntpStrip;
}
  
const NtpMCTruth* MadBase::GetTruth(Int_t i){
  if(!LoadTruth(i)) return 0;
  else return ntpTruth;
}

const NtpMCStdHep* MadBase::GetStdHep(Int_t i){
  if(!LoadStdHep(i)) return 0;
  else return ntpStdHep;
}

const NtpMCTruth* MadBase::GetTruthForReco(Int_t i,Int_t& j){
  if(!LoadTruthForReco(i,j)) return 0;
  else return ntpTruth;
}

const NtpTHSlice* MadBase::GetSliceForRecoTH(Int_t i,Int_t& j){
  if(!LoadSliceForRecoTH(i,j)) return 0;
  else return ntpTHSlice;
}
const NtpMCTruth* MadBase::GetTruthForRecoTH(Int_t i,Int_t& j){
  if(!LoadTruthForRecoTH(i,j)) return 0;
  else return ntpTruth;
}

const NtpTHTrack* MadBase::GetTHTrack(Int_t i){
  if(!LoadTHTrack(i)) return 0;
  else return ntpTHTrack;
}

const NtpTHEvent* MadBase::GetTHEvent(Int_t i){
  if(!LoadTHEvent(i)) return 0;
  else return ntpTHEvent;
}

const NtpTHShower* MadBase::GetTHShower(Int_t i){
  if(!LoadTHShower(i)) return 0;
  else return ntpTHShower;
}

const NtpTHSlice* MadBase::GetTHSlice(Int_t i){
  if(!LoadTHSlice(i)) return 0;
  else return ntpTHSlice;
}

const NtpEMShower* MadBase::GetEMShower(Int_t i){
  if(!LoadEMShower(i)) return 0;
  else return ntpEMShower;
}


#endif // #ifdef madbase_cxx

---