TruthHelper.cxx

Go to the documentation of this file.

#include "DataUtil/TruthHelper.h"

#include "TClonesArray.h"
#include "TArrayI.h"
#include "TArrayF.h"

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"

#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandFitTrackHandle.h"
#include "RecoBase/CandTrackHandle.h"

#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandEventHandle.h"
#include "RecoBase/CandEventListHandle.h"

#include "CandDigit/CandDigitHandle.h"

#include "DataUtil/Truthifier.h"
#include "Digitization/DigiScintHit.h"
#include "Digitization/DigiSignal.h"
#include "Record/SimSnarlRecord.h"
#include "MCNtuple/NtpMCStdHep.h"

ClassImp(TruthHelper)
CVSID("$Id: TruthHelper.cxx,v 1.16 2008/11/05 17:42:27 rodriges Exp $");

TruthHelper::TruthHelper()
{
}

TruthHelper::TruthHelper(const MomNavigator * mom)
  :fMom(mom)
{
}

TruthHelper::~TruthHelper()
{
}


Bool_t TruthHelper::IsNeutrino(const TParticle * part)
{
  if(abs(part->GetPdgCode())==12 ||
     abs(part->GetPdgCode())==14 ||
     abs(part->GetPdgCode())==16)return true;
  return false;
}

Bool_t TruthHelper::IsMuon(const TParticle * part)
{
  if(abs(part->GetPdgCode())==13)return true;
  return false;
}

Bool_t TruthHelper::IsPion(const TParticle * part)
{
  if(abs(part->GetPdgCode())==211)return true;
  return false;
}

Bool_t TruthHelper::IsP0(const TParticle * part)
{
  if(abs(part->GetPdgCode())==111)return true;
  return false;
}

Bool_t TruthHelper::IsGamma(const TParticle * part)
{
  if(abs(part->GetPdgCode())==22)return true;
  return false;
}

Bool_t TruthHelper::IsProton(const TParticle * part)
{
  if(abs(part->GetPdgCode())==2212)return true;
  return false;
}

Bool_t TruthHelper::IsNeutron(const TParticle * part)
{
  if(abs(part->GetPdgCode())==2112)return true;
  return false;
}

Bool_t TruthHelper::IsElectron(const TParticle * part)
{
  if(abs(part->GetPdgCode())==11)return true;
  return false;
}

Bool_t TruthHelper::IsPrimaryShowerPart( const TParticle * part)
{
  if(IsNeutrino(part) || IsMuon(part)) return false;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t status = part->GetStatusCode();
    const TParticle * chainpart = part;
    while (status == 205 || status == 212){
      chainpart =  dynamic_cast<const TParticle*>((*ctca)[chainpart->GetFirstMother()]);
      status = chainpart->GetStatusCode();
    }
    if(status !=1  || IsNeutrino(chainpart) || IsMuon(chainpart)) return false;
  }
  return true;
}

Bool_t TruthHelper::IsDocStatus(const TParticle * part)
{
  // Is the status code "Documentation" or "null"?
  // This should indicate that it was an initial particle
  // into the event generator.
  Int_t status = part->GetStatusCode();
  if( status == 3 || status == 0 ) return true;
  return false;
}


Int_t TruthHelper::GetNeuId(Int_t trackId)
{
  Int_t neuId=-1;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    for( Int_t ind=trackId; ind>-1; ind--){
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)){  //changed from 3 to 0
        neuId=ind;
        return neuId;
      }
    }
  }
  return -1;
}

Int_t TruthHelper::GetNextNeuId(Int_t trackId)
{
  Int_t neuId=-1;
  Int_t siz=-1;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    siz = ctca->GetEntriesFast();
    for( Int_t ind=trackId+1; ind<siz; ind++){
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)) {
        neuId=ind;
        return neuId;
      }
    }
  }
  return siz;
}

Int_t TruthHelper::GetStripNeuIndex(CandStripHandle *csh)
{
   //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  Int_t tracky = -1;
  Bool_t found = false;
  if (ssr){ 
    TIter digitItr(csh->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())) {    
      CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              tracky = abs(hit->TrackId());
              found=true;
            }
          }
        }
      }   
    }
  }
  return tracky;
}

Int_t TruthHelper::GetNeuKinIndex(Int_t trackId)
{
  Int_t neuKinIndex=-1;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    if(trackId<siz && trackId>=0){
      for (Int_t ind=0; ind <= trackId; ++ind) {
        TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
        if(IsNeutrino(part) && IsDocStatus(part)){
          neuKinIndex++;
        } 
      }
    }
  }
  return neuKinIndex;
}


Int_t TruthHelper::GetBestNeuMatch(const CandDigitHandle &cdh)
{
  Int_t bestNeu=-1;
  Double_t maxDE=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  const DigiSignal * signal = truth.GetSignal(cdh);
  if(signal){
    for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
      const DigiScintHit * hit = signal->GetHit(i);
      if(hit)
        {
          if(hit->DE()>maxDE){
            Int_t track = abs(hit->TrackId());
            bestNeu=GetNeuKinIndex(track);
            maxDE=hit->DE();
          }
        }
    }
  }
  return bestNeu;
}

Int_t TruthHelper::GetBestNeuMatch(const CandStripHandle &csh)
{
  Int_t bestNeu=-1;
  Double_t maxDE=0; 
 const Truthifier & truth = Truthifier::Instance(fMom);
 TIter digitItr(csh.GetDaughterIterator());
  while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr()))
    {
      CandDigitHandle tcdh = *cdh;
      const DigiSignal * signal = truth.GetSignal(tcdh);
      if(signal){
        for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
          const DigiScintHit * hit = signal->GetHit(i);
          if(hit)
            {
              if(hit->DE()>maxDE){
                Int_t track = abs(hit->TrackId());
                bestNeu=GetNeuKinIndex(track);
                maxDE=hit->DE();
              }
            }
        }
      }
    }
  return bestNeu;
}
 
Int_t  TruthHelper::GetBestTrackIdMatch(CandTrackHandle& cth)
{
 const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t bestTrack=-1;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI nTrackHits(siz);
    nTrackHits.Reset(0);
    TIter stripitr(cth.GetDaughterIterator());
    while (CandStripHandle* csh = dynamic_cast<CandStripHandle*>(stripitr())) {
      TIter digitItr(csh->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        CandDigitHandle tcdh = *cdh;
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track<siz){
                nTrackHits[track]++;
              }
              else{
                MSG("TruthHelper",Msg::kWarning) << "StdHep size or less than zero" << endl;
              }
              
            }
          }
        }
      }
    }
    Int_t nMax=0;
    for (Int_t ind=0; ind < siz; ++ind) {
      if(nTrackHits[ind]>nMax){
        nMax=nTrackHits[ind];
        bestTrack=ind;
      }
    }   
  }
  return bestTrack;
}
int  TruthHelper::GetBestEventNeuMatch(CandEventHandle &ceh)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t BestNeu=-1;
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr) { 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)) {
        neuIndBuf[nNeutrinos]=ind;                       
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    for(Int_t ind=0; ind< nNeutrinos; ++ind) {
      TIter stripitr(ceh.GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          Bool_t found=false;
          CandDigitHandle tcdh = *cdh;
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId());
                if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]){
                  found=true;
                }
              }
            }
          }     
          if(found) {
            neuEBuf[ind]+=cdh->GetCharge();
          }
        }
      }
    }

    Float_t PEMax=0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind]>PEMax){
        PEMax=neuEBuf[ind];
        BestNeu=neuIndBuf[ind];
      }
    }   
  }
  return BestNeu;
}

Float_t TruthHelper::EventPurity( CandEventHandle& ceh, Int_t neuId)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t hitPE=0;
   
  TIter stripitr(ceh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      CandDigitHandle tcdh = *cdh;
      const DigiSignal * signal = truth.GetSignal(tcdh);
      if(signal){
        for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
          const DigiScintHit * hit = signal->GetHit(i);
          if(hit){
            Int_t track = abs(hit->TrackId());
            if(track>=neuId && track<=nextneuId){
              found=true;
            }
          }
        }
      } 
      if(found) {
        hitPE+=cdh->GetCharge();
        totalPE+= cdh->GetCharge();
      } else {
        totalPE+= cdh->GetCharge();
      }
    }
  }

  if(totalPE>0){
    return hitPE/totalPE; 
  }
  return 0.;
}

Float_t TruthHelper::EventCompleteness( CandEventHandle& ceh, CandSliceHandle &csh)
{
  return EventCompletenessImp(ceh, &csh, false);
}


Float_t TruthHelper::EventCompleteness( CandEventHandle& ceh)
{
  return EventCompletenessImp(ceh, 0, false);
}

Float_t TruthHelper::EventCompletenessAllStrips( CandEventHandle& ceh, 
                                                 CandSliceHandle &csh)
{
  return EventCompletenessImp(ceh, &csh, true);
}


Float_t TruthHelper::EventCompletenessAllStrips( CandEventHandle& ceh)
{
  return EventCompletenessImp(ceh, 0, true);
}

Float_t TruthHelper::EventCompletenessImp( CandEventHandle& ceh, 
                                           CandSliceHandle* csh,
                                           bool useAllStrips)
{
  Float_t truePE=0;
  Float_t evtPE=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId= GetBestEventNeuMatch(ceh);
  Int_t nextneuId= GetNextNeuId(neuId);

  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (!ssr)return 0.0;
  //  const TClonesArray* ctca =
  //   dynamic_cast<const TClonesArray*>
  //   (ssr->FindComponent("TClonesArray","StdHep"));
  
  TIter evtstripitr(ceh.GetDaughterIterator());

  // (PAR) In the case where we're only using strips above the PE cut,
  // we need the minimum strip PE used for the shower reco to get the
  // event completeness, but adding it to CandEvent is too hard, so
  // just use the PE cut used for the shower, if there is one. If
  // there's not, use zero
  CandShowerHandle* cshwh=ceh.GetPrimaryShower();
  Double_t minStripPE= cshwh ? cshwh->GetMinStripPE() : 0.;

  // (PAR) TIter::TIter() is protected, so we can't create a blank
  // iterator that way. This works though
  TIter striplistitr((TIterator*)0);

  if(csh) {
    // (PAR) Only loop through strips in the slice
    striplistitr=csh->GetDaughterIterator();
  } else {
    // (PAR) Loop through all strips in the snarl
    CandRecord *crec = dynamic_cast<CandRecord *>
      (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
    
    CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
      (crec->FindCandHandle("CandStripListHandle"));
    striplistitr=cstriplh->GetDaughterIterator();
  }

  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    // (PAR)
    // Only use slice strips with charge greater than the minimum
    // required for showers, if requested. This means that the event completeness
    // is defined as:
    //
    // (PE in reco event) / (PE in true event ending up in strips above PE cut)
    //
    if ((!useAllStrips) && 
        cstriph->GetCharge(CalDigitType::kPE) < minStripPE) continue;

    TIter digitItr(cstriph->GetDaughterIterator());
    Bool_t found=false;
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              //              TParticle* part = dynamic_cast<TParticle*>((*ctca)[track]);
              if(track> neuId && track<nextneuId){
                //              if(part->GetStatusCode()==1) {
                  found=true;
                  //    }
              }
            }
          }         
        }
      }
    }
    if(found) {
      truePE+=cstriph->GetCharge();
      evtstripitr.Reset();
      Int_t plane=cstriph->GetPlane();
      Int_t strip=cstriph->GetStrip();
      while (CandStripHandle* cevtstriph = dynamic_cast<CandStripHandle*>(evtstripitr())) {
        if(cevtstriph->GetPlane()==plane && cevtstriph->GetStrip()==strip){
          evtPE+=cstriph->GetCharge();
          break;
        }
      }
    }
  }
  if(truePE>0) return (evtPE/truePE);
  return 0.;
}

Int_t TruthHelper::GetBestSliceNeuMatch(CandSliceHandle &csh)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t BestNeu=-1;
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr) { 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)) {
        neuIndBuf[nNeutrinos]=ind;                       
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    for(Int_t ind=0; ind< nNeutrinos; ++ind) {
      TIter stripitr(csh.GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          Bool_t found=false;
          CandDigitHandle tcdh = *cdh;
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId());
                if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]){
                  found=true;
                }
              }
            }
          }     
          if(found) {
            neuEBuf[ind]+=cdh->GetCharge();
          }
        }
      }
    }

    Float_t PEMax=0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind]>PEMax){
        PEMax=neuEBuf[ind];
        BestNeu=neuIndBuf[ind];
      }
    }   
  }
  return BestNeu;
}


Float_t TruthHelper::secondNEU(CandSliceHandle &csh)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t BestNeu=-1;
  Float_t PEMax2 = 0.0;
  Float_t TotPE = 0.0;
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)){  
        neuIndBuf[nNeutrinos]=ind;                       
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      TIter stripitr(csh.GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          Bool_t found=false;
          CandDigitHandle tcdh = *cdh;
          if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found) {
            const DigiSignal * signal = truth.GetSignal(tcdh);
            if(signal){
              for (UInt_t i=0;i<signal->GetNumberOfHits();i++) {
                const DigiScintHit * hit = signal->GetHit(i);
                if(hit){
                  Int_t track = abs(hit->TrackId());
                  if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]) {
                    found=true;
                  }
                }
              }
            }
          }     
          if(found) {
            neuEBuf[ind]+=cdh->GetCharge();
          }
        }  
      }
    }

    Float_t PEMax=0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind]>PEMax){
        PEMax=neuEBuf[ind];
        BestNeu=neuIndBuf[ind];
      }
    }   
    PEMax2=0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind] > PEMax2 && neuEBuf[ind] < PEMax) {
        PEMax2 = neuEBuf[ind];
      }
    }
    for(Int_t ind=0; ind < nNeutrinos; ++ind) {
      TotPE += neuEBuf[ind];
    }
  }
  if(TotPE>0) {
    return (PEMax2/TotPE);
  } else {
    return 0;
  }
}

Int_t TruthHelper::NumNeu(CandSliceHandle &csh)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  //Int_t BestNeu=-1;
  Int_t counter = 0;
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)){  
        neuIndBuf[nNeutrinos]=ind;                       
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      TIter stripitr(csh.GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          Bool_t found=false; 
          CandDigitHandle tcdh = *cdh;
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId()); 
                if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]){
                  found=true;
                }
              }
            }
          }     
          if(found) {
            neuEBuf[ind]+=cdh->GetCharge();
          }
        }
      }
    }
    counter = 0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind]>0){
        counter++;
      }
    }   
  }
  return counter;
}

Float_t TruthHelper::SliceCompleteness(CandSliceHandle & csh)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId=GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t slicePE=0;
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  TIter striplistitr(cstriplh->GetDaughterIterator());
  striplistitr.Reset();
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId){
                found=true;
              }
            }
          }
        }
      }
      if(found) {
        totalPE+= cdh->GetCharge();
      }
    }
  }
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId){
                found=true;
              }
            }
          }
        }
      }
      if(found) {
        slicePE+= cdh->GetCharge();
      }
    }
  }

  if(totalPE>0) return  (slicePE/totalPE);
  return 0.;
}

Int_t TruthHelper::SliceTrueStrip(CandSliceHandle &csh)
{
  //HUGE NOTE:  Since I'm trying to get the number of strips this function is probably correct to 1st order to 99%
  //The reason I say probably is because Strips can be 'polluted' by another neutrino so sometimes a strip that 'truely 
  //belongs to a slice' might have a strip that does belong, but it is polluted by another neutrino -> Neglible effect.
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId = GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Int_t truenumber = 0;
  //First have to find true number of strips so we can make the array.
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  TIter stripitr(cstriplh->GetDaughterIterator());
  stripitr.Reset();
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    Bool_t found=false;
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())) {
      CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found) { 
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++) {
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId) {
                found=true;
              }
            }
          }
        }       
      }
    }    
    if(found) {
      truenumber++;
    }
  }
  return truenumber;
}

Float_t TruthHelper::SlicePurity(CandSliceHandle & csh)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId = GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t hitPE=0;
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found) {
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++) {
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit) {
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId){
                found=true;
              }
            }
          }
        }       
      }
      if(found) {
        hitPE+=cdh->GetCharge();
        totalPE+= cdh->GetCharge();
      } else {
        if(!truth.DigitIsOnlyCrosstalk(tcdh))
          totalPE+= cdh->GetCharge();
      }
    }
  }
  if(totalPE>0){
    return (hitPE/totalPE); 
  }
  return 0;
}

Float_t TruthHelper::SlicePurity_MaxTimeGap(CandSliceHandle & csh, Int_t mincharge)
{
  //Basically a way to measure the purity of a slice after the INITIAL step of slicing, 
  //before adding spectrometer hits and low ph hits.
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId = GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t hitPE=0;
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    if(cstriph->GetPlane() < 121 && cstriph->GetCharge()>=mincharge) {
      TIter digitItr(cstriph->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        Bool_t found=false;
        CandDigitHandle tcdh = *cdh;
        if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit) {
                Int_t track = abs(hit->TrackId());
                if(track>=neuId && track<=nextneuId){
                  found=true;
                }
              }
            }
          }     
        }
        if(found) {
          hitPE+=cdh->GetCharge();
          totalPE+= cdh->GetCharge();
        } else {
          if(!truth.DigitIsOnlyCrosstalk(tcdh))
            totalPE+= cdh->GetCharge();
        }
      }
    }
  }
  if(totalPE>0){
    return (hitPE/totalPE); 
  }
  return 0;
}

Float_t TruthHelper::SliceCompleteness_MaxTimeGap(CandSliceHandle & csh, Int_t mincharge)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId=GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t slicePE=0;
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  TIter striplistitr(cstriplh->GetDaughterIterator());
  striplistitr.Reset();
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    if(cstriph->GetPlane()<121 && cstriph->GetCharge()>=mincharge) {
      TIter digitItr(cstriph->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        Bool_t found=false;
        const CandDigitHandle tcdh = *cdh;
        if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId());
                if(track>=neuId && track<=nextneuId){
                  found=true;
                }
              }
            }
          }
        }
        if(found) {
          totalPE+= cdh->GetCharge();
        }
      }
    }
  }
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    if(cstriph->GetPlane()<121 && cstriph->GetCharge() >=mincharge) {
      TIter digitItr(cstriph->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        Bool_t found=false;
        const CandDigitHandle tcdh = *cdh;
        if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId());
                if(track>=neuId && track<=nextneuId){
                  found=true;
                }
              }
            }
          }
        }
        if(found) {
          slicePE+= cdh->GetCharge();
        }
      }
    }
  }

  if(totalPE>0) return  (slicePE/totalPE);
  return 0.;
}

Int_t TruthHelper::TruthSliceNum() {
  //This function is an 'attempt' to define what the true number of slices should be
  //Obviously this is general and in the end only comparisons of energy spectrums/global things etc,
  //will be the ultimate test of slicing : I rarely use this now
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    TArrayF  neuEBuf2(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)){ 
        neuIndBuf[nNeutrinos]=ind;
        neuEBuf[nNeutrinos]=0;
        neuEBuf2[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    CandRecord *crec = dynamic_cast<CandRecord *>
      (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
    CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
      (crec->FindCandHandle("CandStripListHandle"));
    
    for(Int_t ind=0; ind < nNeutrinos; ++ind) {
      TIter striplistitr(cstriplh->GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
        Bool_t found=false;
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          const CandDigitHandle tcdh = *cdh;
          if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
            const DigiSignal * signal = truth.GetSignal(tcdh);
            if(signal){
              for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
                const DigiScintHit * hit = signal->GetHit(i);
                if(hit && hit->Plane()<121) { //Events should NOT be taken if they interacted wholely in the spectro
                  Int_t track = abs(hit->TrackId());
                  if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]){
                    neuEBuf[ind]++;
                  }
                }
              }
            }
          }
        }
        for(Int_t ind=0; ind < nNeutrinos; ++ind) {
          if(neuEBuf[ind]>0) {
            neuEBuf2[ind]++;
            neuEBuf[ind] = 0;
          } else {
            neuEBuf[ind] = 0;
          }
        }
      }
    }
    
    int counter = 0;
    for(int ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf2[ind] >= 2) counter++;
    }
    
    if(counter > 0) {
      return counter;
    }
    return 0; 
  }
  return 0;
}

Float_t TruthHelper::SliceCompleteness_xtalk(CandSliceHandle & csh)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId=GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t slicePE=0;
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  TIter striplistitr(cstriplh->GetDaughterIterator());
  striplistitr.Reset();
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      const CandDigitHandle tcdh = *cdh;
      if(!found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId) {
                found=true;
              }
            }
          }
        }
      }
      if(found) {
        totalPE+= cdh->GetCharge();
      }
    }
  }
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      const CandDigitHandle tcdh = *cdh;
      if(!found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId) {
                found=true;
              }
            }
          }
        }
      }
      if(found) {
        slicePE+= cdh->GetCharge();
      }
    }
  }

  if(totalPE>0) return (slicePE/totalPE);
  return 0.;
}

Float_t TruthHelper::SlicePurity_xtalk(CandSliceHandle & csh)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId = GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t hitPE=0;
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      CandDigitHandle tcdh = *cdh;
      if(!found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit) {
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId){
                found=true;
              }
            }
          }
        }       
      }
      if(found) {
        hitPE+=cdh->GetCharge();
        totalPE+= cdh->GetCharge();
      } else {
        totalPE+= cdh->GetCharge();
      }
    }
  }
  if(totalPE>0){
    return (hitPE/totalPE); 
  }
  return 0;
}

Double_t TruthHelper::TrueNeuE(CandSliceHandle & csh) 
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  //const Truthifier & truth = Truthifier::Instance(fMom);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  Int_t ind = GetBestSliceNeuMatch(csh);
  if (ssr && ind>=0){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
    if(part->Energy() > 0) {
      MSG("TruthHelper", Msg::kInfo) << "Neu " << ind << " E is: " << part->Energy() << endl;
      return part->Energy();
    } else {
      return 0.0;
    }
  } else {
    return 0.0;
  }
}

Int_t TruthHelper::SliceTrueStripxtalk(CandSliceHandle &csh)
{
  //Same as SliceTrueStrip except now including crosstalk
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId = GetBestSliceNeuMatch(csh);
  Int_t nextneuId = GetNextNeuId(neuId);
  Int_t truenumber1 = 0;
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  TIter stripitr(cstriplh->GetDaughterIterator());
  stripitr.Reset();
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    Bool_t found=false;
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())) {
      CandDigitHandle tcdh = *cdh;
      if(!found) {      
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++) {
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track>=neuId && track<=nextneuId) {
                found=true;
              }
            }
          }
        }       
      }
    }    
    if(found) {
      truenumber1++;
    }
  }
  return truenumber1;
}


Float_t TruthHelper::TrackPurity( CandTrackHandle& cth, Int_t trackId)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Float_t total=0;
  Float_t hits=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    TIter stripitr(cth.GetDaughterIterator());
    while (CandStripHandle* csh = dynamic_cast<CandStripHandle*>(stripitr())) {
      TIter digitItr(csh->GetDaughterIterator());
      Bool_t found=false;
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        CandDigitHandle tcdh = *cdh;
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(track==trackId)found=true;
            }
          }
        }
      } 
      total++;
      if(found) hits++;
    }
  }
  if(total>0){
    return hits/total;
  }
  return 0.0; 
}

Float_t TruthHelper::GetTrackMaxE( CandTrackHandle& cth)
{
  Float_t maxE=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t trackId=GetBestTrackIdMatch(cth);
  Int_t neuId=GetNeuId(trackId);
  Int_t nextneuId = GetNextNeuId(neuId);

  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 


    TIter stripitr(cth.GetDaughterIterator());
    while (CandStripHandle* csh = dynamic_cast<CandStripHandle*>(stripitr())) {
      TIter digitItr(csh->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        CandDigitHandle tcdh = *cdh;
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(abs(hit->ParticleId())==13){
                if(track>neuId && track<nextneuId){
                  if(hit->ParticleEnergy()>maxE)maxE=hit->ParticleEnergy();
                }
              }
            }
          }
        }       
      }
    }
  }
  return maxE; 
}

Float_t TruthHelper::GetTrackMaxE2( CandTrackHandle& cth)
{
  Float_t maxE=0;
  Int_t trackId=GetBestTrackIdMatch(cth);
  Int_t neuId=GetNeuId(trackId);
  Int_t nextneuId = GetNextNeuId(neuId);
  
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){     
    // Get the scint hit array.
    const TObjArray * ScintHitArray = 
      dynamic_cast<const TObjArray*>(ssr->FindComponent(0,"DigiScintHits"));
    if(ScintHitArray==0)return maxE; 
    // Build the array of hits.
    TIter hitarrayIter(ScintHitArray);
    TObject* tobj;
    while( (tobj = hitarrayIter.Next()) ) {
      const DigiScintHit* hit = dynamic_cast<DigiScintHit*>(tobj);
      if(hit){
        Int_t track = abs(hit->TrackId());      
        // look only at muon hits not in veto shield.
        if(abs(hit->ParticleId())==13 && hit->Plane()<500){
          if(track>neuId && track<nextneuId){
            if(hit->ParticleEnergy()>maxE)maxE=hit->ParticleEnergy();
          }
        }
      }
    }
  }
  return maxE; 
}

Float_t TruthHelper::GetTrackMinE2( CandTrackHandle& cth)
{
  Float_t minE=1000;
  Int_t trackId=GetBestTrackIdMatch(cth);
  Int_t neuId=GetNeuId(trackId);
  Int_t nextneuId = GetNextNeuId(neuId);
  
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    
    // Get the scint hit array.
    const TObjArray * ScintHitArray = 
      dynamic_cast<const TObjArray*>(ssr->FindComponent(0,"DigiScintHits"));
    if(ScintHitArray==0)return minE; 
    // Build the array of hits.
    TIter hitarrayIter(ScintHitArray);
    TObject* tobj;
    while( (tobj = hitarrayIter.Next()) ) {
      const DigiScintHit* hit = dynamic_cast<DigiScintHit*>(tobj);
      if(hit){
        Int_t track = abs(hit->TrackId());
        
        // look only at muon hits not in veto shield.
        
        if(abs(hit->ParticleId())==13 && hit->Plane()<500){
          if(track>neuId && track<nextneuId){
            if(hit->ParticleEnergy()<minE)minE=hit->ParticleEnergy();
          }
        }
      }
    }
  }
  return minE; 
}

Float_t TruthHelper::GetTrackMinE( CandTrackHandle& cth)
{
  Float_t minE=10000;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t trackId= GetBestTrackIdMatch(cth);
  Int_t neuId=GetNeuId(trackId);
  Int_t nextneuId = GetNextNeuId(neuId);

  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    
    TIter stripitr(cth.GetDaughterIterator());
    while (CandStripHandle* csh = dynamic_cast<CandStripHandle*>(stripitr())) {
      TIter digitItr(csh->GetDaughterIterator());
      while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
        CandDigitHandle tcdh = *cdh;
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              if(abs(hit->ParticleId())==13){
                if(track>neuId && track<nextneuId){
                  if(hit->ParticleEnergy()<minE)minE=hit->ParticleEnergy();
                }
              }
            }
          }
        }       
      }
    }
  }
  return minE; 
}


Float_t TruthHelper::TrackCompleteness( CandTrackHandle& cth, CandSliceHandle &csh)
{
  Float_t ntruehits=0;
  Float_t ntrkhits=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t trackId=GetBestTrackIdMatch(cth);
  
  TIter trkstripitr(cth.GetDaughterIterator());
  
  TIter striplistitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    Bool_t found=false;
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = hit->TrackId();
              if(track==trackId){
                found=true;
              }
            }
          }         
        }
      }
    }
    if(found) {
      ntruehits++;
      trkstripitr.Reset();
      Int_t plane=cstriph->GetPlane();
      Int_t strip=cstriph->GetStrip();
      while (CandStripHandle* ctrkstriph = dynamic_cast<CandStripHandle*>(trkstripitr())) {
        if(ctrkstriph->GetPlane()==plane && ctrkstriph->GetStrip()==strip){
          ntrkhits++;
          break;
        }
      }
    }
  }
  if(ntruehits>0) return (ntrkhits/ntruehits);
  return 0.;
}

Float_t TruthHelper::TrackCompleteness( CandTrackHandle& cth)
{
  Float_t ntruehits=0;
  Float_t ntrkhits=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t trackId= GetBestTrackIdMatch(cth);
  
  TIter trkstripitr(cth.GetDaughterIterator());
  
  CandRecord *crec = dynamic_cast<CandRecord *>
    (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
//rwh  CandStripSRListHandle *cstriplh = dynamic_cast<CandStripSRListHandle *>
//rwh    (crec->FindCandHandle("CandStripSRListHandle"));
  CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));

  TIter striplistitr(cstriplh->GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    Bool_t found=false;
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = hit->TrackId();
              if(track==trackId){
                found=true;
              }
            }
          }         
        }
      }
    }
    if(found) {
      ntruehits++;
      trkstripitr.Reset();
      Int_t plane=cstriph->GetPlane();
      Int_t strip=cstriph->GetStrip();
      while (CandStripHandle* ctrkstriph = dynamic_cast<CandStripHandle*>(trkstripitr())) {
        if(ctrkstriph->GetPlane()==plane && ctrkstriph->GetStrip()==strip){
          ntrkhits++;
          break;
        }
      }
    }
  }
  if(ntruehits>0) return ntrkhits/ntruehits;
  return 0.;
}

Int_t  TruthHelper::GetBestShowerNeuMatch(CandShowerHandle& csh)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t BestNeu=-1;
  Int_t nNeutrinos=0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)) { 
        neuIndBuf[nNeutrinos]=ind;
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      TIter stripitr(csh.GetDaughterIterator());
      while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
        TIter digitItr(cstriph->GetDaughterIterator());
        while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
          Bool_t found=false;
          CandDigitHandle tcdh = *cdh;
          const DigiSignal * signal = truth.GetSignal(tcdh);
          if(signal){
            for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
              const DigiScintHit * hit = signal->GetHit(i);
              if(hit){
                Int_t track = abs(hit->TrackId());
                if(track>=neuIndBuf[ind] && track <= neuIndBuf[ind+1]){
                  found=true;
                }
              }
            }
          }     
          if(found) {
            neuEBuf[ind]+=cdh->GetCharge();
          }
        }
      }
    }
    Float_t PEMax=0;
    for (Int_t ind=0; ind < nNeutrinos; ++ind) {
      if(neuEBuf[ind]>PEMax){
        PEMax=neuEBuf[ind];
        BestNeu=neuIndBuf[ind];
      }
    }   
  }
  return BestNeu;
}

Float_t TruthHelper::ShowerPurity( CandShowerHandle& csh, Int_t neuId)
{
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t nextneuId = GetNextNeuId(neuId);
  Float_t totalPE=0;
  Float_t hitPE=0;
   
  TIter stripitr(csh.GetDaughterIterator());
  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(stripitr())) {
    TIter digitItr(cstriph->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      Bool_t found=false;
      CandDigitHandle tcdh = *cdh;
      const DigiSignal * signal = truth.GetSignal(tcdh);
      if(signal){
        for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
          const DigiScintHit * hit = signal->GetHit(i);
          if(hit){
            Int_t track = abs(hit->TrackId());
            if(track>=neuId && track<=nextneuId){
              found=true;
            }
          }
        }
      } 
      if(found) {
        hitPE+=cdh->GetCharge();
        totalPE+= cdh->GetCharge();
      } else {
        totalPE+= cdh->GetCharge();
      }
    }
  }

  if(totalPE>0){
    return hitPE/totalPE; 
  }
  return 0.;
}
 
Float_t TruthHelper::ShowerCompleteness( CandShowerHandle& cshwh, CandSliceHandle &csh)
{
  return ShowerCompletenessImp(cshwh, &csh, false);
}

//This has a status code requirement of 1, don't know if that's right or not.
Float_t TruthHelper::ShowerCompleteness( CandShowerHandle& cshwh)
{
  return ShowerCompletenessImp(cshwh, 0, false);
}

Float_t TruthHelper::ShowerCompletenessAllStrips(CandShowerHandle& cshwh, 
                                                 CandSliceHandle &csh)
{
  return ShowerCompletenessImp(cshwh, &csh, true);
}

//This has a status code requirement of 1, don't know if that's right or not.
Float_t TruthHelper::ShowerCompletenessAllStrips( CandShowerHandle& cshwh)
{
  return ShowerCompletenessImp(cshwh, 0, true);
}

Float_t TruthHelper::ShowerCompletenessImp(CandShowerHandle& cshwh, 
                                           CandSliceHandle* csh,
                                           bool useAllStrips)
{
  Float_t truePE=0;
  Float_t shwPE=0;
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t neuId= GetBestShowerNeuMatch(cshwh);
  Int_t nextneuId= GetNextNeuId(neuId);

  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (!ssr)return 0.0;
  const TClonesArray* ctca =
    dynamic_cast<const TClonesArray*>
    (ssr->FindComponent("TClonesArray","StdHep"));
  
  TIter shwstripitr(cshwh.GetDaughterIterator());
  
  // (PAR) TIter::TIter() is protected, so we can't create a blank
  // iterator that way. This works though
  TIter striplistitr((TIterator*)0);

  if(csh) {
    // (PAR) Only loop through strips in the slice
    striplistitr=csh->GetDaughterIterator();
  } else {
    // (PAR) Loop through all strips in the snarl
    CandRecord *crec = dynamic_cast<CandRecord *>
      (fMom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
    
    CandStripListHandle *cstriplh = dynamic_cast<CandStripListHandle *>
      (crec->FindCandHandle("CandStripListHandle"));
    striplistitr=cstriplh->GetDaughterIterator();
  }

  while (CandStripHandle* cstriph = dynamic_cast<CandStripHandle*>(striplistitr())) {
    // (PAR)
    // If requested, only use snarl strips with charge greater than the minimum
    // required for showers. This means that the shower completeness
    // is defined as:
    //
    // (PE in reco shower) / (PE in true shower ending up in strips above PE cut)
    //
    if ((!useAllStrips) && 
        cstriph->GetCharge(CalDigitType::kPE)<cshwh.GetMinStripPE()) continue;
    
    TIter digitItr(cstriph->GetDaughterIterator());
    Bool_t found=false;
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){
      const CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit){
              Int_t track = abs(hit->TrackId());
              TParticle* part = dynamic_cast<TParticle*>((*ctca)[track]);

              if(track>neuId && track<nextneuId && IsPrimaryShowerPart(part)){
                found=true;      
              }
            }
          }
        }
      }
    }
    if(found) {
      truePE+=cstriph->GetCharge();
      shwstripitr.Reset();
      Int_t plane=cstriph->GetPlane();
      Int_t strip=cstriph->GetStrip();
      while (CandStripHandle* cshwstriph = dynamic_cast<CandStripHandle*>(shwstripitr())) {
        if(cshwstriph->GetPlane()==plane && cshwstriph->GetStrip()==strip){
          shwPE+=cstriph->GetCharge();
          break;
        }         
      } 
    }
  }
  if(truePE>0) return shwPE/truePE;
  return 0.;

}

Int_t TruthHelper::Stripxtalk(CandStripHandle *csh) 
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  Int_t flag = 0;
  if (ssr){ 
    TIter digitItr(csh->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())) {    
      CandDigitHandle tcdh = *cdh;
      if(truth.DigitIsOnlyCrosstalk(tcdh)) {
        flag = 1;
      } else {
        flag = 0;
      }
    }
  }
  return flag;
}

Float_t TruthHelper::StripPurity(CandStripHandle *csh)
{
  //Tom Osiecki - osiecki@mail.hep.utexas.edu
  const Truthifier & truth = Truthifier::Instance(fMom);
  Int_t nNeutrinos=0;
  Int_t totalnum = 0;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    TArrayI  neuIndBuf(siz+1);
    TArrayF  neuEBuf(siz);
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && part->GetStatusCode()==0){  
        neuIndBuf[nNeutrinos]=ind;
        neuEBuf[nNeutrinos]=0;
        nNeutrinos++;
      }
    }
    neuIndBuf[nNeutrinos]=siz;
    
    TIter digitItr(csh->GetDaughterIterator());
    while (CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(digitItr())){    
      Bool_t found=false;
      CandDigitHandle tcdh = *cdh;
      if(!truth.DigitIsOnlyCrosstalk(tcdh) && !found){
        const DigiSignal * signal = truth.GetSignal(tcdh);
        if(signal){
          for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
            const DigiScintHit * hit = signal->GetHit(i);
            if(hit) {
              Int_t track = abs(hit->TrackId());
              for(Int_t y=0; y<nNeutrinos; y++) {
                if(track>=neuIndBuf[y] && track<=neuIndBuf[y+1]) {      
                  neuEBuf[y]++;
                }
              }
            }
          }
        }       
      }
    }
    
    for(Int_t y=0; y<nNeutrinos; y++) {
      if(neuEBuf[y]>0) totalnum++;
    }
    
    if(totalnum>0) {
      return totalnum;
    }
  }
  return 0;
}


Int_t TruthHelper::GetClosestNeuVtx(CandRecoHandle& crh)
{
  Int_t bestNeu=-1;
  Float_t vtxU=crh.GetVtxU();
  Float_t vtxV=crh.GetVtxV();
  Float_t vtxX=0.707*(vtxU-vtxV);
  Float_t vtxY=0.707*(vtxU+vtxV);
  Float_t vtxZ=crh.GetVtxZ();
  Float_t closestDist=1e6;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    Int_t siz = ctca->GetEntriesFast();
    for (Int_t ind=0; ind < siz; ++ind) {
      TParticle* part = dynamic_cast<TParticle*>((*ctca)[ind]);
      if(IsNeutrino(part) && IsDocStatus(part)){  //changed from 3 to 0.
        Float_t dist= sqrt((vtxX-part->Vx())*(vtxX-part->Vx())+
                           (vtxY-part->Vy())*(vtxY-part->Vy())+
                           (vtxZ-part->Vz())*(vtxZ-part->Vz()));
        if(dist<closestDist){
          closestDist=dist;
          bestNeu=ind;
        }
      }
    }
  }  
  return bestNeu;
}

const Float_t * TruthHelper::GetP4Neu(Int_t neuId){
 
  for(Int_t i=0;i<4;i++) m_p4Neu[i]=0.;
 if(neuId<0)return m_p4Neu;
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* part = dynamic_cast<TParticle*>((*ctca)[neuId]);
    if(IsNeutrino(part) ){
      m_p4Neu[0]=part->Px();
      m_p4Neu[1]=part->Py();
      m_p4Neu[2]=part->Pz();
      m_p4Neu[3]=part->Energy();
    }
    else{
      
      MSG("TruthHelper",Msg::kWarning) << "Non-neutrino stdhep entry interrigated to obtain Neu 4 vector" << endl;
    } 
  } 
  return m_p4Neu;
}

const Float_t * TruthHelper::GetP4Shw(Int_t neuId){
  for(Int_t i=0;i<4;i++) m_p4Shw[i]=0.;
  if(neuId<0)return m_p4Shw;
  Int_t nextneu = GetNextNeuId(neuId);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* neu = dynamic_cast<TParticle*>((*ctca)[neuId]);
    if(IsNeutrino(neu)){
      for (Int_t i=neuId;i<nextneu ; i++){
        TParticle* part = dynamic_cast<TParticle*>((*ctca)[i]);
        if( !IsNeutrino(part) && !IsMuon(part) && part->GetStatusCode()==1){
          m_p4Shw[0]+=part->Px();
          m_p4Shw[1]+=part->Py();
          m_p4Shw[2]+=part->Pz();
          m_p4Shw[3]+=part->Energy();
          if((part->GetPdgCode()/1000)%10 || part->GetPdgCode()>99999) m_p4Shw[3]-=part->GetMass();     
        }
      }
    }
  }
  return m_p4Shw;
}

const Float_t * TruthHelper::GetP4Mu1(Int_t neuId){

 
  for(Int_t i=0;i<4;i++) m_p4Mu1[i]=0;
 if(neuId<0)return m_p4Mu1;
  Int_t nextneu = GetNextNeuId(neuId);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* neu = dynamic_cast<TParticle*>((*ctca)[neuId]);
    if(IsNeutrino(neu)){
      for (Int_t i=neuId;i<nextneu ; i++){
        TParticle* part = dynamic_cast<TParticle*>((*ctca)[i]);
        if(IsMuon(part)){
          Float_t p[4];
          p[0]=part->Px();
          p[1]=part->Py();
          p[2]=part->Pz();
          p[3]=part->Energy(); 
          if(p[3]>m_p4Mu1[3]){
            m_p4Mu1[0]=p[0];
            m_p4Mu1[1]=p[1];
            m_p4Mu1[2]=p[2];
            m_p4Mu1[3]=p[3];
          }       
        }
      }
    }
  }
  return m_p4Mu1;
}

const Float_t * TruthHelper::GetP4Mu2(Int_t neuId){

  const Float_t * mu1E = GetP4Mu1(neuId);  
  for(Int_t i=0;i<4;i++) m_p4Mu2[i]=0;
 if(neuId<0)return m_p4Mu2;
  Int_t nextneu = GetNextNeuId(neuId);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* neu = dynamic_cast<TParticle*>((*ctca)[neuId]);
    if(IsNeutrino(neu)){
      for (Int_t i=neuId;i<nextneu ; i++){
        TParticle* part = dynamic_cast<TParticle*>((*ctca)[i]);
        if(IsMuon(part) && part->Energy()!=mu1E[3]){
          Float_t p[4];
          p[0]=part->Px();
          p[1]=part->Py();
          p[2]=part->Pz();
          p[3]=part->Energy(); 
          if(p[3]>m_p4Mu2[3]){
            m_p4Mu2[0]=p[0];
            m_p4Mu2[1]=p[1];
            m_p4Mu2[2]=p[2];
            m_p4Mu2[3]=p[3];
          }       
        }
      }
    }
  }
  return m_p4Mu2;
}


const Float_t * TruthHelper::GetP4El1(Int_t neuId){
 
 for(Int_t i=0;i<4;i++) m_p4El1[i]=0; 
 if(neuId<0) return m_p4El1;
  Int_t nextneu = GetNextNeuId(neuId);
  SimSnarlRecord *ssr = 
    dynamic_cast<SimSnarlRecord*>(fMom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca =
      dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    TParticle* neu = dynamic_cast<TParticle*>((*ctca)[neuId]);
    if(IsNeutrino(neu)){
      for (Int_t i=neuId;i<nextneu ; i++){
        TParticle* part = dynamic_cast<TParticle*>((*ctca)[i]);
        if(IsElectron(part)){
          Float_t p[4];
          p[0]=part->Px();
          p[1]=part->Py();
          p[2]=part->Pz();
          p[3]=part->Energy(); 
          if(p[3]>m_p4El1[3]){
            m_p4El1[0]=p[0];
            m_p4El1[1]=p[1];
            m_p4El1[2]=p[2];
            m_p4El1[3]=p[3];
          }       
        }
      }
    }
  }
  return m_p4El1;
}

---