CDTruthifier.cxx

Go to the documentation of this file.

// Program name: CDTruthifier.cxx
//
// Package: CalDetTracker
//
// Purpose: A simple class to hold the MC truth about the hits
//
// Contact: Chris Smith, Ryan Nichol, Leo Jenner or Jeff Hartnell

#include "TClonesArray.h"
#include "TParticle.h"

#include "CalDetPID/CalDetParticleType.h"
#include "MinosObjectMap/MomNavigator.h"
#include "DataUtil/Truthifier.h"

#include "CalDetTracker/CDPIDInfo.h"
#include "CalDetTracker/CDTruthifier.h"
#include "CalDetTracker/CDTruthHitInfo.h"
#include "MessageService/MsgService.h"

using namespace std;

ClassImp(CDTruthifier)

CVSID("$Id: CDTruthifier.cxx,v 1.18 2005/02/02 12:15:14 tagg Exp $");

//......................................................................

CDTruthifier::CDTruthifier()
{
  this->InitialiseVariables();
}

//......................................................................

CDTruthifier::~CDTruthifier()
{
  //nothing
}

//......................................................................

void CDTruthifier::InitialiseVariables()
{
  //nothing
}

//......................................................................

void CDTruthifier::GetTruthHitInfo(const MomNavigator *mom,
                                   TClonesArray* tca)
{ 
  TClonesArray &truthInfo=*tca;   //clonesarray to be filled
  Int_t caCounter=0;
  
  //get the relations between all the data and truth information for MC
  Truthifier truth=Truthifier::Instance(mom);

  //make the complete strip list
  this->MakeCompleteStripList(mom);
  //protect against empty lists
  if (fCompleteStripList.empty()){
    MSG("CDTruthifier",Msg::kWarning)
      <<"No strips in the CompleteStripList, can't get TruthHitInfo"
      <<endl;
    return;
  }

  //Standard pdg codes
  //Int_t electron=11;
  Int_t muon=13;

  //loop over the complete strip list (dsh and/or csh)
  for (map<Int_t,Int_t>::iterator sl=fCompleteStripList.begin();
       sl!=fCompleteStripList.end();++sl){

    //get the key, plane and strip
    Int_t key=sl->first;
    //Build18BitPlnStripKey=(GetPlane()<<bitsPlexIdStrip) | GetStrip();
    Int_t bitsPlexIdStrip=8;
    Int_t mask=static_cast<Int_t>(pow(2.,bitsPlexIdStrip)-1);
    Int_t strip=mask & key;
    Int_t plane=key>>bitsPlexIdStrip;

    //create a new CDTruthHitInfo in the TClonesArray
    new(truthInfo[caCounter]) CDTruthHitInfo(plane,strip);

    MSG("CDTruthifier",Msg::kVerbose)
      <<"pl="<<plane<<", st="<<strip<<endl;

    Double_t totalEnergyDep=0;
    Int_t mainParticle=0;
    Double_t mainPathLength=0;
    Double_t mainParticleEnergy=-1;
    Int_t digiScintHitCounter=0;
    Double_t mainT1=0;
    Double_t mainX1=0;
    Double_t mainY1=0;
    Double_t mainZ1=0;
    Double_t mainT2=0;
    Double_t mainX2=0;
    Double_t mainY2=0;
    Double_t mainZ2=0;
    Double_t initialEarliestT1=1e50;
    Double_t earliestT1=initialEarliestT1;
    Double_t latestT2=-1e-9;
    Bool_t foundFirstMainParticle=false;
    Bool_t foundFirstParticle=false;

    //loop over the digiScintHits
    for (multimap<Int_t,DigiScintHit*>::iterator hitIter=
           fAllDigiScintHits.equal_range(key).first;
         hitIter!=fAllDigiScintHits.equal_range(key).second;
         ++hitIter){

      //count the digiScintHits
      digiScintHitCounter++;

      //get the pointer to the dsh
      DigiScintHit* dsh=hitIter->second;

      //sanity check the digiScintHits actually belong to the strip
      //and are not just connected by shared pmts or something
      if (!(plane==dsh->Plane() &&
            strip==dsh->Strip())) {
        MSG("CDTruthifier",Msg::kWarning)
          <<"Ahhhhhhhh, wrong plane and strip"<<endl;
      }
      
      //get the dsh info
      const Double_t energyDep=dsh->DE();
      const Double_t pathLength=dsh->DS();
      const Double_t energyParticle=dsh->ParticleEnergy();
      const Int_t particleId=dsh->ParticleId();
      const Int_t trackId=dsh->TrackId();
      const Double_t T1=dsh->T1();
      const Double_t T2=dsh->T2();
      
      //sum up the energy depositions of the different particles
      //in this strip
      totalEnergyDep+=energyDep;

      //look for particles other than electrons and muons
      if (particleId!=11 && particleId!=-11 && 
          particleId!=13 && particleId!=-13 &&
          particleId!=211 && particleId!=-211 && //pi 
          particleId!=2212 && particleId!=-2212 && //proton
          particleId!=1002001000 && //what is this???
          particleId!=1003001000 && //what is this???
          particleId!=1004001000 && //what is this???
          particleId!=1004002000){ //what is this??? alpha nucleus?
        MSG("CDTruthifier",Msg::kInfo)
          <<endl<<endl
          <<"Found particles other than expected"<<endl
          <<digiScintHitCounter
          <<" DE="<<energyDep*1000
          <<", DS="<<pathLength*100
          <<", pId="<<particleId
          <<", pEn="<<energyParticle*1000
          <<", tEn="<<totalEnergyDep*1000
          <<", tId="<<trackId
          <<", mP="<<mainParticle
          <<", mPL="<<mainPathLength*100<<endl<<endl;
      }

      //set the times
      //check if earliest or latest
      if (T1<earliestT1) earliestT1=T1;
      if (T2>latestT2) latestT2=T2;

      //this section sets the variables in the case where there
      //is no main particle
      //they are set to the first particle in the list, dodgy???        
      if (!foundFirstParticle){
        foundFirstParticle=true;

        //set main particle and energy
        mainParticle=particleId;
        mainParticleEnergy=energyParticle;
        
        //set the ?1s
        mainT1=dsh->T1();
        mainX1=dsh->X1();
        mainY1=dsh->Y1();
        mainZ1=dsh->Z1();
        
        //set the path length
        mainPathLength=pathLength;
        //set the ?2s
        mainT2=dsh->T2();
        mainX2=dsh->X2();
        mainY2=dsh->Y2();
        mainZ2=dsh->Z2();
      }

      //check if the "main" particle
      //need a better way to do this! Need to know the mother particle
      if (abs(particleId)==muon){

        static Int_t firstPlane=-1;
        static Int_t firstStrip=-1;

        //for a given particle you can get multiple DigiScintHits 
        //per strip so need to be smart in how this is handled
        if (!foundFirstMainParticle){
          foundFirstMainParticle=true;
 
          //only set main particle and energy once
          mainParticle=particleId;
          mainParticleEnergy=energyParticle;

          //only set the ?1s once too
          mainT1=dsh->T1();
          mainX1=dsh->X1();
          mainY1=dsh->Y1();
          mainZ1=dsh->Z1();

          //set the path length
          //in the case of multiple hits this will be summed below
          mainPathLength=pathLength;
          //set the ?2s
          mainT2=dsh->T2();
          mainX2=dsh->X2();
          mainY2=dsh->Y2();
          mainZ2=dsh->Z2();

          //set control variables
          firstPlane=plane;
          firstStrip=strip;
        }
        else{

          //set the message level variable
          Msg::LogLevel_t msgLevel=Msg::kDebug;

          //check that the plane and strip are the same as the
          //first muon
          if (!(plane==firstPlane && strip==firstStrip)){
            //set the message level so the next message prints
            msgLevel=Msg::kWarning;
            MSG("CDTruthifier",Msg::kWarning)
              <<endl
              <<"Plane and strip are not the same as the first main"
              <<"particle, First (pl,st)=("<<firstPlane<<","<<firstStrip
              <<") New=("<<plane<<","<<strip<<")"
              <<endl;
          }

          MSG("CDTruthifier",msgLevel)
            <<"Values already set:"<<endl
            <<"DigiScintHit="<<digiScintHitCounter
            <<endl
            <<"  Current:"
            <<" (pl,st)=("<<plane<<","<<strip<<")"
            <<" PL="<<mainPathLength*100
            <<", PE="<<mainParticleEnergy*1000
            <<endl
            <<"  T1="<<mainT1*1e9
            <<", X1="<<mainX1*100
            <<", Y1="<<mainY1*100
            <<", Z1="<<mainZ1*100
            <<endl
            <<"  T2="<<mainT2*1e9
            <<", X2="<<mainX2*100
            <<", Y2="<<mainY2*100
            <<", Z2="<<mainZ2*100
            <<endl
            <<"  New:"
            <<" (pl,st)=("<<dsh->Plane()<<","<<dsh->Strip()<<")"
            <<" PL="<<pathLength*100
            <<", PE="<<energyParticle*1000
            <<", tEnDep="<<totalEnergyDep*1000<<"("<<energyDep*1000<<")"
            <<", partId="<<particleId
            <<", trkId="<<trackId
            <<endl
            <<"  T1="<<dsh->T1()*1e9
            <<"  X1="<<dsh->X1()*100
            <<", Y1="<<dsh->Y1()*100
            <<", Z1="<<dsh->Z1()*100
            <<endl
            <<"  T2="<<dsh->T2()*1e9
            <<"  X2="<<dsh->X2()*100
            <<", Y2="<<dsh->Y2()*100
            <<", Z2="<<dsh->Z2()*100
            <<endl;

          //only change the values of those below if on same strip
          //and plane as the first main particle to be found
          if (plane==firstPlane && strip==firstStrip){
            //this is not an "else" of the above because of the 
            //large message

            MSG("CDTruthifier",Msg::kDebug)
              <<"Comparing path length for multiple DigiScintHits"<<endl
              <<" for a single particle in a single strip:"<<endl
              <<"DigiScintHit="<<digiScintHitCounter
              <<endl
              <<"  Current:"
              <<" (pl,st)=("<<plane<<","<<strip<<")"
              <<" PL="<<mainPathLength*100
              <<endl
              <<"  New:"
              <<" (pl,st)=("<<dsh->Plane()<<","<<dsh->Strip()<<")"
              <<" PL="<<(mainPathLength+pathLength)*100
              <<"("<<pathLength*100<<")"
              <<endl;

            //set the path length
            //in the case of multiple hits it is summed
            mainPathLength+=pathLength;
            
            //the ?2s should be set to those of the last dsh  
            //assuming ordered in time
            if (mainT2<=dsh->T2()){
              mainT2=dsh->T2();
              mainX2=dsh->X2();
              mainY2=dsh->Y2();
              mainZ2=dsh->Z2();
            }
            else{
              MSG("CDTruthifier",Msg::kWarning)
                <<endl
                <<"(pl,st)=("<<dsh->Plane()<<","<<dsh->Strip()<<")"
                <<" Multiple dsh for particle in given strip"<<endl
                <<" are not time ordered, this may give false"
                <<" reconstruction. Beware!"<<endl;
            }
          }
        }
      }

      //general message
      MSG("CDTruthifier",Msg::kVerbose)
        <<digiScintHitCounter
        <<" DE="<<energyDep*1000
        <<", DS="<<pathLength*100
        <<", pId="<<particleId
        <<", pEn="<<energyParticle*1000
        <<", tEn="<<totalEnergyDep*1000
        <<", tId="<<trackId
        <<", mP="<<mainParticle
        <<", mPL="<<mainPathLength*100<<endl;
    }

    Int_t finalPmtTruth1=DigiSignal::kUnknown;
    Int_t finalPmtTruth2=DigiSignal::kUnknown;
    Int_t vaChip1=-1;
    Int_t vaChip2=-1;

    //create a pointer for current csh
    CandStripHandle* csh=0;

    //check that csh exists (it often will not)
    if ((csh=fAllCandStripHandles[key])){

      //get iterator over csh daughters
      TIter cshItr(csh->GetDaughterIterator());
      
      //counter for the digits
      Int_t digitCounter=0;

      //now loop over the canddigits
      while(CandDigitHandle* cdh=dynamic_cast<CandDigitHandle*>
            (cshItr.Next())){
        //count the digits
        digitCounter++;
        
        DigiSignal::DigiPmtTruth_t pmtTruth=truth.
          GetSignalTruthFlags(*cdh);
        
        //set the final pmt truth
        //this is hard to deal on a whole strip basis
        //you could have crosstalk on one end but genuine on 
        //the other
        //this just sets all the bits that it finds
        //how to interpret it can be done in the analysis
        if (cdh->GetPlexSEIdAltL().GetEnd()==StripEnd::kEast){
          finalPmtTruth1|=pmtTruth;//OR it with itself
          vaChip1=cdh->GetChannelId().GetVaChip();
        }
        else if (cdh->GetPlexSEIdAltL().GetEnd()==StripEnd::kWest){
          finalPmtTruth2|=pmtTruth;
          vaChip2=cdh->GetChannelId().GetVaChip();
        }
        
        //stop assigning the value of 65535
        if (vaChip2>2) vaChip2=-2;//use -2 to distinguish between when
        if (vaChip1>2) vaChip1=-2;//not set at all and when qie
        
        MSG("CDTruthifier",Msg::kVerbose)
          <<"Dig="<<digitCounter
          <<", va1="<<vaChip1
          <<", va2="<<vaChip2
          <<", pmtTr="<<pmtTruth
          <<", T="<<cdh->GetTime()*1e9
          <<", el="<<ElecType::AsString(cdh->GetChannelId().
                                        GetElecType())
          <<endl;
      }//end of while
      
      //reset all the energies if there is no genuine hit
      //you can get given the energy on a strip that caused a 
      //crosstalk hit, thus double counting
      //NOTE: This must be inside the "if csh" since
      //might have genuine digiScintHit but no digitisation
      //start by taking the OR of the two stripends
      if (!(((finalPmtTruth1 | finalPmtTruth2) & 
             DigiSignal::kGenuine)==DigiSignal::kGenuine)){
        mainParticle=0;
        mainParticleEnergy=-1;
        mainPathLength=0;
        mainT1=0;
        mainX1=0;
        mainY1=0;
        mainZ1=0;
        mainT2=0;
        mainX2=0;
        mainY2=0;
        mainZ2=0;
        totalEnergyDep=0;
      }

    }//end of if csh
    else{
      //set truth as genuine for strips with no digitisation
      finalPmtTruth1=DigiSignal::kGenuine;
      finalPmtTruth2=DigiSignal::kGenuine;

      //leave the vaChips as initialised, no need to do it here
      //vaChip1=-1;
      //vaChip2=-1;
    }

    MSG("CDTruthifier",Msg::kVerbose)
      <<"Final:"
      <<", part="<<mainParticle
      <<", pmtTr1="<<finalPmtTruth1
      <<", pmtTr2="<<finalPmtTruth2
      <<", mainPartEn="<<mainParticleEnergy*1000
      <<", totEn="<<totalEnergyDep*1000
      <<endl;
    
    //change the earliest time to some
    if (earliestT1==initialEarliestT1) earliestT1=-1e-9;

    //get a pointer to the current thi in the tclonesarray
    CDTruthHitInfo* thi=(CDTruthHitInfo*)(truthInfo.At(caCounter));

    //fill the current CDTruthHitInfo 
    thi->SetPmtTruth1(finalPmtTruth1);
    thi->SetPmtTruth2(finalPmtTruth2);
    thi->SetMainParticle(mainParticle);
    thi->SetTotalEnDep(totalEnergyDep);
    thi->SetMainPartEn(mainParticleEnergy);
    thi->SetMainPathLength(mainPathLength);
    thi->SetVaChip1(vaChip1);
    thi->SetVaChip2(vaChip2);
    thi->SetMainT1(mainT1);
    thi->SetMainX1(mainX1);
    thi->SetMainY1(mainY1);
    thi->SetMainZ1(mainZ1);
    thi->SetMainT2(mainT2);
    thi->SetMainX2(mainX2);
    thi->SetMainY2(mainY2);
    thi->SetMainZ2(mainZ2);
    thi->SetEarliestT1(earliestT1);
    thi->SetLatestT2(latestT2);
    thi->SetNumDigiScintHits(digiScintHitCounter);

    //increment the clones array counter
    caCounter++;
  }
  //kUnknown          = 0x00, // Can happen: known cause is 0 charge
  // but random noise passes sparsification thresh.
  //kGenuine          = 0x01, //A real hit from real gosh-darn physics.
  //kDarkNoise        = 0x02, // Dark noise from the tube.
  //kFibreLight       = 0x04, // Random [singles] from a fibre. 
  //kCrosstalk        = 0x08, // A generic crosstalk hit.
  //kCrosstalkOptical = 0x10, // A crosstalk hit, specifically from
//kLeakFromNextBucket = 0x20, // The next bucket leaked charge back to this one
//kLeakFromPrevBucket = 0x40  // The prev bucket leaked charge forward to this one
}


//......................................................................

void CDTruthifier::RecoStdHep(const MomNavigator* mom,
                              CDPIDInfo* pidInfo) const
{
  //get the sim snarl record
  SimSnarlRecord *ssr= 
    dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
  if (ssr){ 
    const TClonesArray* ctca=dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    
    Int_t num=ctca->GetEntries();
    MSG("CDTruthifier",Msg::kVerbose)
      <<"StdHep Entries="<<num<<endl;

    Int_t pdgCode=0;
    Float_t particleMomentum=0;

    for(Int_t i=0;i<num;i++){
      TParticle* part=dynamic_cast<TParticle*>((*ctca)[i]);
      
      MSG("CDTruthifier",Msg::kVerbose)
        <<"i="<<i<<", status code="<<part->GetStatusCode()
        <<", pdg code="<<part->GetPdgCode()
        <<", en="<<part->Energy()<<", P="<<part->P()<<endl;

      //set the variables
      pdgCode=part->GetPdgCode();
      particleMomentum=part->P();

      if (i==1){
        MSG("CDTruthifier",Msg::kWarning)
          <<"More than one particle in simsnarlrecord"<<endl
          <<"i="<<i<<", status code="<<part->GetStatusCode()
          <<", pdg code="<<part->GetPdgCode()
          <<", en="<<part->Energy()<<", P="<<part->P()<<endl;
      }
    }

    //set the pid info
    if (pidInfo){
      pidInfo->SetTrueParticleMomentum(particleMomentum);
      pidInfo->SetPIDType(this->PdgCode2CalDetParticleType(pdgCode));
    }
    else{
      MSG("CDTruthifier",Msg::kWarning)
        <<"CDPIDInfo pointer="<<pidInfo<<", can't set truth info"<<endl;
    }
  }
  else{
    MSG("CDTruthifier",Msg::kWarning)
      <<"CDTruthifier::RecoStdHep:: No simsnarl record found"<<endl;
  }
}

//......................................................................

Float_t CDTruthifier::ParticleMomentum(const MomNavigator* mom) const
{
  //get the sim snarl record
  SimSnarlRecord *ssr= 
    dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));

  Float_t particleMomentum=0;
  
  if (ssr){ 
    const TClonesArray* ctca=dynamic_cast<const TClonesArray*>
      (ssr->FindComponent("TClonesArray","StdHep"));
    
    Int_t num=ctca->GetEntries();
    MSG("CDTruthifier",Msg::kVerbose)
      <<"StdHep Entries="<<num<<endl;

    for(Int_t i=0;i<num;i++){
      TParticle* part=dynamic_cast<TParticle*>((*ctca)[i]);
      
      MSG("CDTruthifier",Msg::kVerbose)
        <<"i="<<i<<", status code="<<part->GetStatusCode()
        <<", pdg code="<<part->GetPdgCode()
        <<", en="<<part->Energy()<<", P="<<part->P()<<endl;

      //set the variables
      particleMomentum=part->P();

      if (i==1){
        MSG("CDTruthifier",Msg::kWarning)
          <<"More than one particle in simsnarlrecord"<<endl
          <<"i="<<i<<", status code="<<part->GetStatusCode()
          <<", pdg code="<<part->GetPdgCode()
          <<", en="<<part->Energy()<<", P="<<part->P()<<endl;
      }
    }
  }
  else{
    MSG("CDTruthifier",Msg::kWarning)
      <<"No simsnarl record found"<<endl;
  }
  
  return particleMomentum;
}

//......................................................................

Int_t CDTruthifier::PdgCode2CalDetParticleType(Int_t pdgCode) const
{
  Int_t cdParticleType=CalDetParticleType::kUnknown;

  if (abs(pdgCode)==13){
    cdParticleType=CalDetParticleType::kMuon;
  }
  else if (abs(pdgCode)==211){
    cdParticleType=CalDetParticleType::kPion;
  }
  else if (abs(pdgCode)==11){
    cdParticleType=CalDetParticleType::kElectron;
  }
  else if (abs(pdgCode)==2211){
    cdParticleType=CalDetParticleType::kProton;
  }
  else if (abs(pdgCode)==321){
    cdParticleType=CalDetParticleType::kKaon;
  }
  else{
    cdParticleType=CalDetParticleType::kUnknown;
    MSG("CDTruthifier",Msg::kWarning)
      <<"Particle PDG code not recognised="<<pdgCode<<endl;
  }
    
  return cdParticleType;
}

//......................................................................

void CDTruthifier::ExtractAllDigiScintHits(const MomNavigator *mom)
{
  //clear the multimap
  fAllDigiScintHits.clear();

  SimSnarlRecord* simsnarl=0;
  TObject* tobj=0;
  TIter fragiter=mom->FragmentIter();
  
  // Get the simsnarl.
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }
  
  // Verify there IS a simsnarl.
  if(!simsnarl) {
    MSG("CDTruthifier",Msg::kError)<<"No SimSnarl found"<<endl;
    return;
  }

  // Get the DigiScintHits.
  const TObjArray* hitList=dynamic_cast<const TObjArray*>
    (simsnarl->FindComponent(0,"DigiScintHits"));
  if(hitList==0) {
    MSG("CDTruthifier",Msg::kError) << "Can't find scint hit array.\n"; 
  }

  MSG("CDTruthifier",Msg::kVerbose)
    <<"Found digiScintHits, size="<<hitList->GetSize()<<endl;

  for (Int_t i=0;i<hitList->GetSize();i++){
    DigiScintHit* hit=0;
    if ((hit=dynamic_cast<DigiScintHit*>(hitList->At(i)))){
      
      MSG("CDTruthifier",Msg::kVerbose)
        <<"i="<<i<<", plane="<<hit->Plane()<<endl;

      //get the key
      Int_t key=hit->StripEndId().Build18BitPlnStripKey();

      fAllDigiScintHits.insert(pair<Int_t,DigiScintHit*>(key,hit));
    }
  }  

  MSG("CDTruthifier",Msg::kVerbose)
    <<"mm size="<<fAllDigiScintHits.size()<<endl;

  Int_t counter=0;
  for (multimap<Int_t,DigiScintHit*>::iterator digiScintHit=
         fAllDigiScintHits.equal_range(40).first;
       digiScintHit!=fAllDigiScintHits.equal_range(40).second;
       ++digiScintHit){
    counter++;
    MSG("CDTruthifier",Msg::kInfo)
      <<"i="<<counter<<", plane="<<digiScintHit->second->Plane()
      <<", DE="<<digiScintHit->second->DE()
      <<", PID="<<digiScintHit->second->ParticleId()<<endl;
  }
}

//......................................................................

void CDTruthifier::ExtractAllCandStripHandles(const MomNavigator *mom)
{
  //clear the map
  fAllCandStripHandles.clear();

  CandRecord* candRecord=0;
  TObject* tobj=0;
  TIter fragiter=mom->FragmentIter();

  // Get the candRecord.
  while( ( tobj = fragiter.Next() ) ) {
    candRecord = dynamic_cast<CandRecord*>(tobj);
    if(candRecord) break;
  }

  // Verify there IS a candRecord.
  if(!candRecord) {
    MSG("CDTruthifier",Msg::kError)<<"No CandRecord found"<<endl;
    return;
  }

  //Get CandStripListHandle
  CandStripListHandle *cslh = dynamic_cast<CandStripListHandle*>
    (candRecord->FindCandHandle("CandStripListHandle"));

  //verify there is a cslh
  if(cslh==0) {
    MSG("CDTruthifier", Msg::kFatal) 
      << "No CandStripListHandle in MOM" << endl;
    return;
  }

  //Loop over all candstrips in list
  TIter stripItr(cslh->GetDaughterIterator());
  //TObject* obj;
  while(CandStripHandle* csh=dynamic_cast<CandStripHandle*>
        (stripItr.Next())){ 
    //get the key
    Int_t key=csh->GetStripEndId().Build18BitPlnStripKey();

    //add the csh to the map
    fAllCandStripHandles.insert(pair<Int_t,CandStripHandle*>(key,csh));
  }

  MSG("CDTruthifier",Msg::kVerbose)
    <<"cshs size="<<fAllCandStripHandles.size()<<endl;

}

//......................................................................

void CDTruthifier::MakeCompleteStripList(const MomNavigator *mom)
{
  //clear the list from previous events
  fCompleteStripList.clear();
  
  //extract the digiScintHits and csh from mom
  this->ExtractAllDigiScintHits(mom);
  this->ExtractAllCandStripHandles(mom);
  
  //loop over all the digiScintHits
  //could try and speed this up by jumping to the end of each
  //key section using upper_bound(key)
  Int_t lastKey=-1;
  for (multimap<Int_t,DigiScintHit*>::iterator digiScintHit=
         fAllDigiScintHits.begin();
       digiScintHit!=fAllDigiScintHits.end();
       ++digiScintHit){

    //get the key
    Int_t key=digiScintHit->first;

    //check if key has already been added to map
    //doing this is quicker than searching the map
    //everytime
    if (static_cast<Int_t>(key)==lastKey) continue;
    lastKey=static_cast<Int_t>(key);
    
    //add it to the list or just increment the number of times
    //it was accessed
    fCompleteStripList[key]++;
  }

  //loop over all the csh
  for (map<Int_t,CandStripHandle*>::iterator cshIter=
         fAllCandStripHandles.begin();
       cshIter!=fAllCandStripHandles.end();
       ++cshIter){

    //get the key
    Int_t key=cshIter->first;

    //add it to the list or just increment the number of times
    //it was accessed
    fCompleteStripList[key]++;
  }

  MSG("CDTruthifier",Msg::kDebug)
    <<"Complete list size="<<fCompleteStripList.size()
    <<", all csh size="<<fAllCandStripHandles.size()
    <<", all dsh size="<<fAllDigiScintHits.size()<<endl;
}

//......................................................................

---