Truthifier.cxx

Go to the documentation of this file.

#include "Truthifier.h"
#include "MessageService/MsgService.h"
#include "RecoBase/CandStripHandle.h"
#include <TClonesArray.h>

CVSID( " $Id: Truthifier.cxx,v 1.14 2006/05/19 14:55:42 rhatcher Exp $");
ClassImp(Truthifier)





// ----------------------------------------------------------------------- //
const Truthifier& Truthifier::Instance( const MomNavigator* mom, Bool_t supressErrors )
{

  // Deal with the case of a null mom.
  static Truthifier nullTruth(0);
  if(mom==0) return nullTruth;

  // Look for the simSnarl in Mom.
  SimSnarlRecord* simsnarl = 0;
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }
  
  if(simsnarl==0) return nullTruth; // No good.

  // look for the Truthifier in the simSnarl.
  const TObject* o;
  o = simsnarl->FindTemporary("Truthifier","Truthifier");
  if(o) {
    // Found one.
    Truthifier* ct = const_cast<Truthifier*>(dynamic_cast<const Truthifier*>(o));
    if(ct){
      if(ct->IsValid()) return *ct;
      else {
        ct->Reset(mom,supressErrors);
        return *ct;
        // it's not valid. Tell the record to drop it.
        // whooops, this is illegal.
        //simsnarl->RemoveTemporary(o);
      }
    }    
  }

  // Ok, now there's no truthifier.
  // Make one.
  Truthifier* t = new Truthifier(mom, supressErrors);

  // Save it.
  simsnarl->AdoptTemporary(t);
  
  // return it.
  return *t;
}

// ----------------------------------------------------------------------- //
void Truthifier::Reset( const MomNavigator* mom, Bool_t supressErrors ) 
{

  fMom = mom;
  fRddb = 0;
  fStdHep = 0;
  fNeutrinoList.clear();
  fTrackIndexToNeutrinoIndex.clear();
  fNeutrinoIndexToTrackIndex.clear();
  fTrackToHit.clear();
  fHitToSignal.clear();
  fSignalToDigit.clear();
  fSignals.clear(); 
  fValid=false;

  if(fMom ==0 ) return;
  
  // Find the simsnarl.
  SimSnarlRecord* simsnarl = 0;
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();

  // Get the simsnarl.
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }
  if(simsnarl==0) { MSG("DataUtil",Msg::kWarning) << "Can't find simsnarl.\n"; return; };
  
  // Get the scint hit array.
  fScintHitArray = 
    dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiScintHits"));
  if(fScintHitArray==0) { 

    if(!supressErrors) MSG("DataUtil",Msg::kWarning) << "Can't find scint hit array.\n"; 

  } else {
  
    // Build the array of hits.
    TIter hitarrayIter(fScintHitArray);
    while( (tobj = hitarrayIter.Next()) ) {
      const DigiScintHit* scinthit = dynamic_cast<DigiScintHit*>(tobj);
      if(scinthit)
        fTrackToHit[scinthit->TrackId()].push_back(scinthit);
    }

  }

  // Get the signal array.
  const TObjArray* signalArray = 
    dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiSignal"));
  if(signalArray==0) { 

    if(!supressErrors) MSG("DataUtil",Msg::kWarning) << "Can't find signal array.\n";
  
  } else {

    // Build the array of signals.
    TIter signalIter(signalArray);
    while( (tobj = signalIter.Next()) ) {
      const DigiSignal* signal = dynamic_cast<DigiSignal*>(tobj);
      if(signal) {
        fSignals[signal->GetId()] = signal;
        for(UInt_t i=0;i<signal->GetNumberOfHits();i++) {
          const DigiScintHit* hit = signal->GetHit(i);
          if(hit) 
            fHitToSignal[hit].push_back(signal);
        }
      }
    }
    if(!supressErrors) MSG("DataUtil",Msg::kDebug) << "Got " << fSignals.size() << " signals." << std::endl;

  }

  // Get the StdHep list.
  fStdHep = dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","StdHep"));
  if(fStdHep ==0) {
   if(!supressErrors)  MSG("DataUtil",Msg::kWarning) << "Can't find StdHep list.\n";
  } else {
    // Look through initial pariticle list. Associate track ID with originating neutrinos
    Int_t siz = fStdHep->GetEntriesFast();
    fTrackIndexToNeutrinoIndex.resize(siz);
    fNeutrinoList.clear();
    Int_t neu = 0;
    for(Int_t ind=0; ind<siz; ind++) {
      TParticle* part = dynamic_cast<TParticle*>((*fStdHep)[ind]);
      
      // Is this a lepton?
      if(abs(part->GetPdgCode())>=12 &&
         abs(part->GetPdgCode())<=16 
         ) {
        // 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 ) {
          // We got a neutrino event!
          neu=ind;
          fNeutrinoList.push_back(ind);
        }
      }
      fTrackIndexToNeutrinoIndex[ind]=neu;
      fNeutrinoIndexToTrackIndex[neu].push_back(ind);
    }    
  }

  // Get the NeuKin list.
  fNeuKins = dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","NeuKinList"));
  
  // Find the RDDB.
  const RawHeader* rawhead =0;
  fRddb =0;
  fragiter.Reset();
  while( ( tobj = fragiter.Next() ) ) {
    RawRecord* rawRec = dynamic_cast<RawRecord*>(tobj);
    if(rawRec) {
      rawhead = rawRec->GetRawHeader();
      TIter itr = rawRec->GetRawBlockIter();
      RawDataBlock* rawBlk;
      while ( ( rawBlk = dynamic_cast<RawDataBlock*>(itr.Next()) ) ) {
        fRddb = dynamic_cast<RawDigitDataBlock*>(rawBlk);
        if(fRddb) {
          // Got one!
          TIter diter = fRddb->GetDatumIter();
          RawDigit* rawDigit;
          diter.Reset();
          while( ( rawDigit = (RawDigit*) diter.Next() ) ) {
            // Relate this digit to the signal
            const DigiSignal* signal = GetSignal(rawDigit);
            if(signal) {
              fSignalToDigit[signal] = rawDigit;
            }
          }
        }
      }
    }
  } // Ugly ugly ugly

  // See if we're good enough.
  fValid = true;


  if(fRddb==0) {
    //if(!supressErrors) MSG("DataUtil",Msg::kError) << "Couldn't find RDDB. Truthifier is invalid. \n";
    fValid = false;
  }
  if(fTrackToHit.size()==0) {
    //if(!supressErrors) MSG("DataUtil",Msg::kError) << "No hits found. Truthifier is invalid. \n";
    fValid = false;
  }
  if(fHitToSignal.size()==0) {
    //if(!supressErrors) MSG("DataUtil",Msg::kError) << "No map hit->signal. Truthifier is invalid. \n";
    fValid = false;
  }
  if(fSignalToDigit.size()==0) {
    //if(!supressErrors) MSG("DataUtil",Msg::kError) << "No map signal->digit found. Truthifier is invalid. \n";
    fValid = false;
  }
  if(fSignals.size()==0) {
    //if(!supressErrors)  MSG("DataUtil",Msg::kError) << "No signals found. Truthifier is invalid. \n";
    fValid = false;
  }  
  if ((!fValid)&&(!supressErrors)) {
    MAXMSG("DataUtil",Msg::kWarning,5) 
      << "Truthifier built but is invalid. Not MC or null snarl." << endl;
  }
}

// Utility.
std::vector<Int_t> Truthifier::GetListOfNeutrinoIndecies() const
{
  return fNeutrinoList;
}

const TParticle*   Truthifier::GetStdHepParticle(Int_t index) const
{
  if(fStdHep) {
    TParticle* part = dynamic_cast<TParticle*>((*fStdHep)[index]);
    return part;
  }
  return 0;
}

Int_t   Truthifier::GetStdHepSize() const
{
  if(fStdHep) {
    return fStdHep->GetEntries();
  }
  return 0;
}

const REROOT_NeuKin* Truthifier::GetNeuKin(Int_t i) const
{
  if(fNeuKins) {
    return dynamic_cast<const REROOT_NeuKin*>((*fNeuKins)[i]);
  }
  return 0;
}

const REROOT_NeuKin* Truthifier::GetNeuKinOfParticle(Int_t index) const
{
  // Find matching neutrino
  for(Int_t i=fNeutrinoList.size()-1;i>=0;i--) {
    if(fNeutrinoList[i] <= index) {
      return GetNeuKin(i);
    }
  }
  return 0;
}

Int_t   Truthifier::GetScintHitListSize() const
{
  if(fScintHitArray) {
    return fScintHitArray->GetEntriesFast();
  }
  return 0;
}

const DigiScintHit* Truthifier::GetScintHit(Int_t index) const
{
  if(index >= GetScintHitListSize()) return 0;
  return dynamic_cast<const DigiScintHit*>( (*fScintHitArray)[index] );
}



const RawDigit*   Truthifier::GetRawDigit(const CandDigitHandle& cdh) const
{
  int digitIndex = cdh.GetRawDigitIndex();
  const RawDigit* digit = 0;
  if(fRddb) digit = fRddb->At(digitIndex);
  return digit;
}


// FORWARD TRACKING
// 
// Routines to go from Particle->Digit
//

Int_t Truthifier::GetNeutrinoOfTrack( Int_t track ) const
{
  return fTrackIndexToNeutrinoIndex[abs(track)];
}


// ----------------------------------------------------------------------- //
const std::vector<const DigiSignal*>& 
Truthifier::GetSignalVector( const DigiScintHit* hit ) const
{
  HitToSignal_t::const_iterator it = fHitToSignal.find(hit);
  if(it!=fHitToSignal.end()) return it->second;

  static std::vector<const DigiSignal*> nullResult;
  return nullResult;
}

// ----------------------------------------------------------------------- //
const RawDigit* 
Truthifier::GetRawDigit( const DigiSignal* signal ) const
{
  SignalToDigit_t::const_iterator it = fSignalToDigit.find(signal);
  if(it!=fSignalToDigit.end()) return it->second;
  return 0;
}

// ----------------------------------------------------------------------- //
const std::vector<const RawDigit*>& 
Truthifier::GetRawDigitVector( const DigiScintHit* hit ) const
{
  static std::vector<const RawDigit*> output;
  output.clear();
  const std::vector<const DigiSignal*>& signalList
    = GetSignalVector(hit);
  for(UInt_t i=0;i<signalList.size();i++) {
    output.push_back(GetRawDigit(signalList[i]));
  }
  return output;
}


// ----------------------------------------------------------------------- //
template<class T>
std::vector<const T*> Truthifier::GetCandidatesFromParticle(int ParticleId, 
                                                              const CandHandle* list) const
{
  std::vector<const T*> output;

  // Loop throuth the list.
  TIter it = list->GetDaughterIterator();
  TObject* obj;
  while( (obj=it.Next()) ) {
    const T* thing = dynamic_cast<const T*>(obj);
    if(thing) {
      const CandHandle* handle = dynamic_cast<const CandHandle*>(thing);
      if(handle)
        if(IsCandidateFromParticle(ParticleId,handle)) output.push_back(thing);
    }
  }; 
  
  return output;
}

// ----------------------------------------------------------------------- //
template<class T>
std::vector<const T*> Truthifier::GetCandidatesFromTrack(int TrackId, 
                                                         const CandHandle* list) const
{
  std::vector<const T*> output;

  // Loop throuth the list.
  TIter it = list->GetDaughterIterator();
  TObject* obj;
  while( (obj=it.Next()) ) {
    const T* thing = dynamic_cast<const T*>(obj);
    if(thing) {
      const CandHandle* handle = dynamic_cast<const CandHandle*>(thing);
      if(handle)
        if(IsCandidateFromTrack(TrackId,handle)) output.push_back(thing);
    }
  }; 
  
  return output;
}


/*
  // For me.
void Truthifier::Test() const 
{
  std::vector<const CandDigitHandle*> output;
  
  if(fMom==0) return;
  const CandHandle* list = DataUtil::GetCandidate<CandHandle>
    (fMom,"CandStripSRListHandle","CandStripSRList");

  CandDigitsFrom(output,list);
  MSG("DataUtil",Msg::kInfo) << "total cand digits:" << output.size() << endl;

  std::vector<const CandStripHandle*> outlist;
  outlist = GetCandidatesFromParticle<CandStripHandle>(13,list);
  
  MSG("DataUtil",Msg::kInfo) << "total muon strips:" << outlist.size() << endl;  
}
*/

// BACKWARD TRACKING
// 
// Routines to go from Digit->Particle
//


// ----------------------------------------------------------------------- //
const DigiSignal* Truthifier::GetSignal( const RawDigit* digit ) const
{

  if(digit==0) return NULL;

  const RawDigiDigitMixIn* ddigit = dynamic_cast<const RawDigiDigitMixIn*>(digit);
  if(ddigit==0) return NULL;

  UInt_t id = ddigit->GetDigiSignalIndex();
  
  SignalMap_t::const_iterator it = fSignals.find(id);
  if(it!=fSignals.end()) return it->second;

  return NULL;
}

// ----------------------------------------------------------------------- //
const DigiSignal* Truthifier::GetSignal( const CandDigitHandle& cdh ) const
{
  int digitIndex = cdh.GetRawDigitIndex();
  const RawDigit* digit = 0;
  if(fRddb) digit = fRddb->At(digitIndex);
  return GetSignal(digit);
}


// ----------------------------------------------------------------------- //
DigiSignal::DigiPmtTruth_t
Truthifier::GetSignalTruthFlags( const CandDigitHandle& cdh ) const
{
  const DigiSignal* signal = GetSignal(cdh);
  if(signal==0) return DigiSignal::kUnknown;
  return signal->GetTruth();
}

// ----------------------------------------------------------------------- //
DigiSignal::DigiPmtTruth_t 
Truthifier::GetSignalTruthFlags( const RawDigit* digit ) const
{
  const DigiSignal* signal = GetSignal(digit);
  if(signal==0) return DigiSignal::kUnknown;
  return signal->GetTruth();
}

// ----------------------------------------------------------------------- //
Bool_t Truthifier::DigitIsOnlyCrosstalk( const CandDigitHandle& cdh ) const
{
  DigiSignal::DigiPmtTruth_t flags = GetSignalTruthFlags(cdh);
  if(  ( (flags & DigiSignal::kCrosstalk) || (flags & DigiSignal::kCrosstalkOptical) )
       && ( ! (flags & DigiSignal::kGenuine) ) ) {
    return true;
  }
  return false;
}
 
// ----------------------------------------------------------------------- //
Bool_t Truthifier::DigitIsOnlyCrosstalk( const RawDigit* digit ) const
{
  DigiSignal::DigiPmtTruth_t flags = GetSignalTruthFlags(digit);
  if(  ( (flags & DigiSignal::kCrosstalk) || (flags & DigiSignal::kCrosstalkOptical) )
       && ( ! (flags & DigiSignal::kGenuine) ) ) {
    return true;
  }
  return false;

}


// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsSignalFromParticle( Int_t particleID,  const DigiSignal* signal ) const
{
  if(! signal) return 0;
  Float_t total = 0;
  Float_t totalMatch = 0;
  UInt_t n = signal->GetNumberOfHits();
  for(UInt_t i=0; i<n; i++) {
    const DigiScintHit* hit = signal->GetHit(i);
    Float_t weight = signal->GetHitWeight(i);
    if(hit->ParticleId() == particleID)  totalMatch += weight;
    total += weight;
  }
  if(total<=0) return 0;  // Guard against FPEs
  return totalMatch/total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsSignalFromTrack( Int_t trackID,  const DigiSignal* signal ) const
{
  if(! signal) return 0;
  Float_t total = 0;
  Float_t totalMatch = 0;
  UInt_t n = signal->GetNumberOfHits();
  for(UInt_t i=0; i<n; i++) {
    const DigiScintHit* hit = signal->GetHit(i);
    Float_t weight = signal->GetHitWeight(i);
    if(hit->TrackId() == trackID)  totalMatch += weight;
    total += weight;
  }
  if(total<=0) return 0;  // Guard against FPEs
  return totalMatch/total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsSignalFromNeutrino( Int_t neutrinoID,  const DigiSignal* signal, Double_t max_t ) const
{
  if(! signal) return 0;
  Float_t total = 0;
  Float_t totalMatch = 0;
  UInt_t n = signal->GetNumberOfHits();
  for(UInt_t i=0; i<n; i++) {
    const DigiScintHit* hit = signal->GetHit(i);
    Float_t weight = signal->GetHitWeight(i);
    total += weight;
    int trackId = abs(hit->TrackId());

    if(fTrackIndexToNeutrinoIndex[trackId]==neutrinoID) {
      if(hit->T1()<=max_t)
        totalMatch += weight;
    }
  }
  if(total<=0) return 0;  // Guard against FPEs
  return totalMatch/total;
}


// ----------------------------------------------------------------------- //
const DigiScintHit* 
Truthifier::GetBiggestHit( const CandDigitHandle& cdh ) const
{
  const DigiSignal* signal = GetSignal(cdh);
  if(!signal) return 0;
  return signal->GetBiggestHit();
}
// ----------------------------------------------------------------------- //
const DigiScintHit* 
Truthifier::GetBiggestHit( const RawDigit* digit ) const
{
  const DigiSignal* signal = GetSignal(digit);
  if(!signal) return 0;
  return signal->GetBiggestHit();
}



// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromParticle( Int_t particleID,  const RawDigit* digit ) const
{
  return IsSignalFromParticle(particleID,GetSignal(digit));
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromParticle( Int_t particleID,  const CandDigitHandle& cdh ) const
{
  return IsSignalFromParticle(particleID,GetSignal(cdh));
}


// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromTrack( Int_t trackID,  const RawDigit* digit ) const
{
  return IsSignalFromTrack(trackID,GetSignal(digit));
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromNeutrino( Int_t neutrinoID,  const RawDigit* digit, Double_t max_t ) const
{
  return IsSignalFromNeutrino(neutrinoID,GetSignal(digit),max_t);
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromTrack( Int_t trackID,  const CandDigitHandle& cdh ) const
{
  return IsSignalFromTrack(trackID,GetSignal(cdh));
}


// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsDigitFromNeutrino( Int_t neutrinoID,  const CandDigitHandle& cdh, Double_t max_t ) const
{
  return IsSignalFromNeutrino(neutrinoID,GetSignal(cdh),max_t);
}



// ----------------------------------------------------------------------- //
PlexStripEndId Truthifier::BestSEIdOfDigit(const CandDigitHandle* cdh) const
{
  // Get list of the stripends.
  std::map<PlexStripEndId,float> ends;
  std::map<PlexStripEndId,float>::iterator endit;

  if(!cdh) {
    MSG("DataUtil",Msg::kError) << "BestSEIdOfDigit was handed null CandDigitHandle*" << endl;
    return PlexStripEndId();
  }

  const DigiSignal* signal = GetSignal(*cdh);
  if(signal) {
    for(UInt_t i=0; i<signal->GetNumberOfHits(); i++) {
      PlexStripEndId seid = signal->GetHit(i)->StripEndId();
      float weight = signal->GetHitWeight(i);
      ends[seid] += weight;
    }
  } else return PlexStripEndId();

  // Find the best one.
  PlexStripEndId best;
  float bestpe = 0;
  for(endit = ends.begin(); endit!= ends.end(); endit++) {
    if(endit->second > bestpe) {
      bestpe = endit->second;
      best = endit->first;
    }
  }
  
  if(best.IsValid()) best.SetEnd(cdh->GetPlexSEIdAltL().GetEnd());
  
  return best;
}

// ----------------------------------------------------------------------- //
void
Truthifier::GetAllScintHits(const RawDigit* digit,       
                            std::set<const DigiScintHit*>& outSet) const
{
  const DigiSignal* signal = GetSignal(digit);
  if(!signal) return;
  UInt_t n = signal->GetNumberOfHits();
  for(UInt_t i=0; i<n; i++) {
    const DigiScintHit* hit = signal->GetHit(i);
    outSet.insert(hit);
  }
}

// ----------------------------------------------------------------------- //
void
Truthifier::GetAllScintHits(const CandDigitHandle& cdh,   
                            std::set<const DigiScintHit*>& outSet) const
{
  const DigiSignal* signal = GetSignal(cdh);
  if(!signal) return;
  UInt_t n = signal->GetNumberOfHits();
  for(UInt_t i=0; i<n; i++) {
    const DigiScintHit* hit = signal->GetHit(i);
    outSet.insert(hit);
  }
}

// ----------------------------------------------------------------------- //
void
Truthifier::GetAllScintHits(const CandHandle* candidate,   
                            std::set<const DigiScintHit*>& outSet) const
{
  if(!candidate) return;
  std::vector<const CandDigitHandle*> digits;
  CandDigitsFrom(digits,candidate);

  for(UInt_t i=0;i<digits.size();i++) 
    GetAllScintHits(*(digits[i]),outSet);
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::GetTotalEnergy( const RawDigit* digit, Int_t particleID ) const 
{
  // Get the total energy of all scinthits that lent energy to this digit.
  // If particleID is non-zero, only add energy from that particle type
  // (-ve or +ve of that type)
 
  std::set<const DigiScintHit*> hits;
  std::set<const DigiScintHit*>::iterator it;
  Float_t total = 0;
  GetAllScintHits(digit,hits);
  for(it = hits.begin(); it!= hits.end(); it++) {
    const DigiScintHit* hit = *it;
    if(hit)
      if((particleID==0) || (abs(hit->ParticleId()) == particleID))
        total += hit->DE();
  }

  return total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::GetTotalEnergy( const CandDigitHandle& cdh, Int_t particleID ) const 
{
  // Get the total energy of all scinthits that lent energy to this digit.
  // If particleID is non-zero, only add energy from that particle type
  // (-ve or +ve of that type)
 
  std::set<const DigiScintHit*> hits;
  std::set<const DigiScintHit*>::iterator it;
  Float_t total = 0;
  GetAllScintHits(cdh,hits);
  for(it = hits.begin(); it!= hits.end(); it++) {
    const DigiScintHit* hit = *it;
    if(hit)
      if((particleID==0) || (abs(hit->ParticleId()) == particleID))
        total += hit->DE();
  }

  return total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::GetTotalEnergy( const CandHandle* candidate, Int_t particleID ) const
{
  // Get the total energy of all scinthits that lent energy to this digit.
  // If particleID is non-zero, only add energy from that particle type
  // (-ve or +ve of that type)
 
  std::set<const DigiScintHit*> hits;
  std::set<const DigiScintHit*>::iterator it;
  Float_t total = 0;
  GetAllScintHits(candidate,hits);
  for(it = hits.begin(); it!= hits.end(); it++) {
    const DigiScintHit* hit = *it;
    if(hit)
      if((particleID==0) || (abs(hit->ParticleId()) == particleID))
        total += hit->DE();
  }

  return total;
}



// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsCandidateFromParticle( Int_t particleID,  const CandHandle* candidate ) const
{
  if(!candidate) return 0;
  std::vector<const CandDigitHandle*> digits;
  CandDigitsFrom(digits,candidate);

  // Test each digit. to see if it satisfies.
  Float_t total = 0;
  for(UInt_t i=0;i<digits.size();i++) 
    total += IsDigitFromParticle(particleID,*digits[i]);
  
  return total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsCandidateFromTrack( Int_t trackID,  const CandHandle* candidate ) const
{
  if(!candidate) return 0;
  std::vector<const CandDigitHandle*> digits;
  CandDigitsFrom(digits,candidate);

  // Test each digit. to see if it satisfies.
  Float_t total = 0;
  for(UInt_t i=0;i<digits.size();i++) 
    total += IsDigitFromTrack(trackID,*digits[i]);
  
  return total;
}

// ----------------------------------------------------------------------- //
Float_t 
Truthifier::IsCandidateGenuine( const CandHandle* candidate ) const
{
  if(!candidate) return 0;
  std::vector<const CandDigitHandle*> digits;
  CandDigitsFrom(digits,candidate);

  // Test each digit. to see if it satisfies.
  Float_t total = 0;
  for(UInt_t i=0;i<digits.size();i++) {
    DigiSignal::DigiPmtTruth_t flags = GetSignalTruthFlags(*digits[i]);       
    if(flags & DigiSignal::kGenuine) {
      total++;
    }
  }

  return total/(float)(digits.size());
}

// ----------------------------------------------------------------------- //
void Truthifier::Print(const Option_t*) const
{

  // Go through all tracks.
  TrackToHit_t::const_iterator hitItr = fTrackToHit.begin();
  for( ; hitItr!=fTrackToHit.end(); hitItr++) {
    std::cout << "+--- TRACK " << hitItr->first << " ---------\n";
    
    // For each scint hit
    for(UInt_t ihit=0; ihit<hitItr->second.size(); ihit++) {
      const DigiScintHit* hit =  (hitItr->second)[ihit];
      std::cout << "|+-- DigiScintHit: " << hit->AsString() << std::endl;
      
      // For each signal
      HitToSignal_t::const_iterator htsItr = fHitToSignal.find(hit);
      if(htsItr != fHitToSignal.end()) {
        for(UInt_t isig = 0; isig<htsItr->second.size(); isig++) {
          
          const DigiSignal* signal = (htsItr->second)[isig];
          SignalToDigit_t::const_iterator stdItr = fSignalToDigit.find(signal);
          if(stdItr != fSignalToDigit.end()) {
            const RawDigit* digit = stdItr->second;
            std::cout << " |-- RawDigit: " << *digit << std::endl;
            std::cout << " |   DigiSignal: " 
                      << signal->GetCharge() / Munits::fC 
               << " fC" 
                      << "  Truth: " << AsString(signal->GetTruth()) << std::endl;      
          }
        }
        
      }
      
    }
  }
}


// Private functions.
// 
//

// ----------------------------------------------------------------------- //
void Truthifier::CandDigitsFrom(std::vector<const CandDigitHandle*>& output, 
                                const CandHandle* list) const
{
  if(!list) return;
 
  const CandDigitHandle* cdh = 
    dynamic_cast<const CandDigitHandle*>(list);
 
  if(cdh){
    output.push_back(cdh);
    return;
  }

  TIter it = list->GetDaughterIterator();
  TObject* obj;
  while( (obj=it.Next()) ) {
    const CandHandle* newHandle = dynamic_cast<const CandHandle*>(obj);
    if(newHandle)
      CandDigitsFrom(output,newHandle);
  }; 
}

---