PTSimApplication.cxx

Go to the documentation of this file.

//
// PTSimApplication
//
// PTSimApplication is a concrete MCAppUser application for use in simulating
// transport of particles through the MINOS detector geometries.
//
// Author:  S. Kasahara 05/04

#include <TVirtualMC.h>
#include <TArrayI.h>
#include <TArrayD.h>
#include <TGeoManager.h>
#include <TGeoNode.h>
#include <TClonesArray.h>
#include <TDatabasePDG.h>

#include "GeoGeometry/GeoStripNode.h"
#include "Digitization/DigiScintHit.h"
#include "Util/UtilPDG.h"
#include "MCApplication/MCApplication.h"
#include "ParticleTransportSim/PTSimApplication.h"
#include "ParticleTransportSim/PTSimStack.h"
#include "G3PTSim/TGeant3/TGeant3.h"

#include "MessageService/MsgService.h"

CVSID("$Id: PTSimApplication.cxx,v 1.67 2009/06/22 03:51:01 kasahara Exp $");

ClassImp(PTSimApplication)

//_____________________________________________________________________________
PTSimApplication::PTSimApplication(const char* name, const char* title) 
  : MCAppUser(name,title),fSnarl(-1),fEvent(-1),fPrimary(-1),fNSnarls(0),
    fNEvents(0),fPTSimStack(0),fHitArray(0),fStdHepArray(0),fMCHit(),
    fExitLiner(false),fLogLevel(kInfo),fLastStackSize(0),fNHitBegEvt(0),
    fNSensitiveMedId(0),fCurrentPath(""),fRkVLevel(0),fRkVMinDistInCM(-1.),
    fRkVFactor(-1.),fRockdEdXMin(-1.) {
  // Normal constructor

  MSG("PTSim",Msg::kDebug) << "PTSimApplication normal ctor @ " 
                           << this << endl;

  fPTSimStack = new PTSimStack();

}

//_____________________________________________________________________________
PTSimApplication::~PTSimApplication() {
  // delete all the owned sub-objects

  MSG("PTSim",Msg::kDebug) << "PTSimApplication dtor @ " 
                           << this << endl;
 
  if ( fPTSimStack ) delete fPTSimStack; fPTSimStack = 0;
  ResetEvtStdHepMap();

}

//_____________________________________________________________________________
void PTSimApplication::GeneratePrimaries() {
  // Fill the user stack with primary particles

  MSG("PTSim",Msg::kDebug) << "PTSimApplication::GeneratePrimaries " 
                           << fSnarl << "/" << fEvent << endl;

  // Clear stack
  fPTSimStack->Reset();
  
  // Push primaries to stack for this event
  Int_t nevent = fEvtStdHepMap.size();
  if ( fEvent >= nevent ) {
    MSG("PTSim",Msg::kError) << "GeneratePrimaries called with Event "
                             << fEvent << " but only " << nevent
                             << " particles in stack. Abort." << endl;
    abort();
  }

  std::vector<TParticle*>& stdheplist = fEvtStdHepMap[fEvent];

  TVector3 pol;
  for ( unsigned int is = 0; is < stdheplist.size(); is++ ) {
    TParticle* prim = stdheplist[is];

    Int_t toBeDone = 0; // only track final states
    if ( prim->GetStatusCode() == 1 ) {
      if ( abs(prim->GetPdgCode()) == 311 ) {
        // If Monte Carlo type is TGeant3, convert K0/K0bar to K0L/K0S 
        // so that Geant3 can handle them
        TGeant3* geant3 = dynamic_cast<TGeant3*>(gMC);
        if ( geant3 ) {
          Double_t randno = MCApplication::Instance().GetRandom() -> Rndm();
          if ( randno <= 0.5 ) prim -> SetPdgCode(130); // K0L
          else prim -> SetPdgCode(310); // K0S
        }
      }
      Int_t primPdg = prim->GetPdgCode();
      if ( primPdg != 28 && primPdg != 0 ) {
        // particles 28 (reggeon) and 0 (rootino) unknown to MC
        if ( gMC->IdFromPDG(primPdg) <= 0 ) {
          // particle unknown to monte carlo and not reggeon or rootino
          // print warning
          const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
          std::string pname = "???";
          if (dbpdg.GetParticle(primPdg)) 
                                pname = dbpdg.GetParticle(primPdg)->GetName();
          MSG("PTSim",Msg::kWarning) << "Particle of type "
                                     << primPdg << "/" << pname.c_str()  
                                     << " not defined to Monte Carlo."
                                     << " Skipped!" << endl;
        }
        else toBeDone = 1;
      }
    }

    Int_t parent = prim->GetMother(0); // parent trkid, -1 => no parent track
    prim->GetPolarisation(pol);
    Int_t trkId = -1; // filled by stack
    // kPPrimary is used for creator process VMC code
    fPTSimStack -> PushTrack(toBeDone,parent,prim->GetPdgCode(),
                             prim->Px(),prim->Py(),prim->Pz(),prim->Energy(),
                             prim->Vx(),prim->Vy(),prim->Vz(),prim->T(),
                             pol.X(),pol.Y(),pol.Z(),kPPrimary,trkId,
                             prim->GetWeight(),prim->GetStatusCode());
     
  }
  
  if ( fLogLevel <= Msg::kDebug ) {
    MSG("PTSim",Msg::kDebug) << "PTSimApplication::GeneratePrimaries "
                             << fSnarl << "/" << fEvent << " pushed "
                             << fPTSimStack->GetNtrack() << " particles w/"
                             << fPTSimStack->GetNprimary() << " primaries "
                             << "to stack." << endl;
    fPTSimStack->Print();
  }
  
  return;

}

//_____________________________________________________________________________
void PTSimApplication::BeginEvent() {    
  // User actions at beginning of event

  fEvent++; // event in snarl
  fNEvents++; // total number of events
  
  MSG("PTSim",Msg::kDebug) << "  *** PTSimApplication::BeginEvent: Snarl " 
                           << fSnarl << " Event " << fEvent 
                           << " ***  " << endl;
  fPrimary = -1;

  bool toswim = false;
  const_cast<MCApplication&>(MCApplication::Instance()).SwitchMedia(toswim);

  // Log the sensitive media here.  This is used to improve the
  // performance during stepping, since the sensitive media is likely
  // to be a small subset of all the media.
  fSensitiveMedId.clear();
  TIter next(gGeoManager->GetListOfMedia());
  TGeoMedium* med = 0;
  fNSensitiveMedId = 0;
  while ((med = (TGeoMedium*)next())) {
    if ( med->GetParam(0) > 0 ) {
      // sensitive medium
      std::string medName = med -> GetName();
      if ( medName.find("_SWIM") == std::string::npos ) {
        fSensitiveMedId.push_back(med->GetId());
        fNSensitiveMedId++;
      }
    }
  }

  fLogLevel = MsgService::Instance()->GetStream("PTSim")->GetLogLevel();
  fLastStackSize = fPTSimStack->GetNtrack();
  fNHitBegEvt    = fHitArray->GetEntriesFast();
  
    
}

//_____________________________________________________________________________
void PTSimApplication::BeginPrimary() {    
  // User actions at beginning of a primary track

  fPrimary++;
  if ( fLogLevel <= Msg::kDebug ) {
    const TParticle* primary = fPTSimStack->GetCurrentTrack();
    Int_t pdg = primary->GetPdgCode();
    std::string pname = primary->GetName();
    Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
 
    MSG("PTSim",Msg::kDebug) << "PTSimApplication::BeginPrimary " 
                             << fSnarl << "/" << fEvent << "/"
                             << fPrimary <<  " trkno " << trkno << " pdg "
                             << pdg << "/" << pname.c_str() <<  endl;
  }  
  
}

//_____________________________________________________________________________
void PTSimApplication::PreTrack() {    
  // User actions at beginning of each track

  fExitLiner = false; // reset exit flag

  if ( fLogLevel <= Msg::kDebug ) {
    Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
    TLorentzVector pos;
    gMC->TrackPosition(pos);
    TLorentzVector mom;
    gMC->TrackMomentum(mom);
    Int_t pdg = gMC->TrackPid();
    const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
    std::string pname = "???";
    if (dbpdg.GetParticle(pdg)) pname = dbpdg.GetParticle(pdg)->GetName();
    Int_t parentno = fPTSimStack->GetCurrentParentTrackNumber();

    if ( fLogLevel == Msg::kDebug ) {
      if ((fPTSimStack -> GetParticle(trkno) -> GetProcess()) 
                                                         == kPPrimary) {
        MSG("PTSim",Msg::kDebug) << "PTSimApplication::PreTrack " 
          << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
          << " parent " << parentno   
          << "\n   " << pdg << "/" << pname.c_str() 
          << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
          << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
          << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
          << " " << mom.E() << endl;
      }
    }
    else {
      MSG("PTSim",Msg::kVerbose) << "PTSimApplication::PreTrack " 
        << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
        << " parent " << parentno   
        << "\n   " << pdg << "/" << pname.c_str() 
        << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
        << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
        << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
        << " " << mom.E() << endl;
    }
  }

}

//_____________________________________________________________________________
void PTSimApplication::Stepping() {    
  // User actions at each step


  if ( gMC -> IsNewTrack() ) {
    Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
    if ( fRkVLevel > 1 || (fRkVLevel == 1  
                           && (fPTSimStack -> GetParticle(trkno) -> GetProcess()) == kPPrimary) ) {
      if ( IsRockVetoed() ) {
        if ( fLogLevel <= Msg::kDebug ) {
          TLorentzVector pos;
          gMC->TrackPosition(pos);
          TLorentzVector mom;
          gMC->TrackMomentum(mom);
          Int_t pdg = gMC->TrackPid();
          const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
          std::string pname = "???";
          if (dbpdg.GetParticle(pdg)) pname=dbpdg.GetParticle(pdg)->GetName();
          Int_t parentno = fPTSimStack->GetCurrentParentTrackNumber();
          MSG("PTSim",Msg::kDebug) 
            << "Track is in rock and below required energy, call StopTrack\n" 
            << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
            << " parent " << parentno   
            << "\n   " << pdg << "/" << pname.c_str() 
            << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
            << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
            << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
            << " " << mom.E() << endl;
        }
        gMC -> StopTrack();
        return;
      }
    }
  }

  if ( fLogLevel <= Msg::kVerbose ) {

    TLorentzVector msgpos;
    gMC->TrackPosition(msgpos);
    TLorentzVector msgmom;
    gMC->TrackMomentum(msgmom);

    MSG("PTSim",Msg::kVerbose) << "Stepping:\n"
    << "particle Id " << gMC->TrackPid() << " in vol " 
    <<  gMC->CurrentVolName() 
    << "\npos " << msgpos.X() << " " << msgpos.Y() << " " << msgpos.Z() << " " 
    << msgpos.T() << "\nmom " << msgmom.X() << " " << msgmom.Y() << " " 
    << msgmom.Z() << " " << msgmom.E() << "\neloss " << gMC->Edep() 
    << ", step " << gMC->TrackStep() << endl;

    TArrayI processlist;
    Int_t nprocess = gMC->StepProcesses(processlist);
    std::string processnamelist = "Step Processes:\n";
    for (int ip = 0; ip < nprocess; ip++ ) {
      processnamelist +=  std::string(TMCProcessName[processlist[ip]]) + "/";
    }
    MSG("PTSim",Msg::kVerbose) << processnamelist.c_str() << endl;
  
    MSG("PTSim",Msg::kVerbose) << "Track Is:" 
          << (gMC->IsTrackEntering()    ? "Entering/" : "" )  
          << (gMC->IsTrackExiting()     ? "Exiting/" : "" )  
          << (gMC->IsTrackOut()         ? "Out/" : "" )  
          << (gMC->IsTrackStop()        ? "Stop/" : "" )  
          << (gMC->IsNewTrack()         ? "New/" : "" )  
          << (gMC->IsTrackInside()      ? "Inside/" : "" )  
          << (gMC->IsTrackDisappeared() ? "Disappeared/" : "" )  
          << (gMC->IsTrackAlive()       ? "Alive/" : "" )  
          << endl;
  }
  
  Int_t nsec = gMC->NSecondaries();
  if ( nsec != 0 ) {
    if ( fLogLevel <= Msg::kVerbose ) {
      std::string secondarystr;
      if ( nsec == 1 ) 
        secondarystr = " secondary generated in this step.";
      else 
        secondarystr = " secondaries generated in this step.";
      MSG("PTSim",Msg::kVerbose) << nsec << secondarystr.c_str() << endl;
    }
  
    if ( fPTSimStack->GetStdHepSaveSibling() ) {
      // Record secondary siblings.
      // if we have recorded a decay, increment fLastStackSize by 1
      // because the decay parent at point of decay will have been
      // pushed to stack
      TMCProcess prodprocess = gMC->ProdProcess(0);
      if ( prodprocess == kPDecay ) fLastStackSize++;
    
      if ( nsec != (fPTSimStack->GetNtrack() - fLastStackSize) ) {
        // Sanity check, difference between total stack size and
        // last step stack size should be equal to number of secondaries
        // + the parent particle stored at the decay point.
        // If not, shut down.
        MSG("PTSim",Msg::kError) 
        << "Number of secondaries generated in this step " << nsec 
        << "\nnot equal to difference in stack size since last addition."
        << "Was: " << fLastStackSize << " Now: " << fPTSimStack->GetNtrack()
        << " Abort." << endl;
        abort();
      }
      for ( Int_t isec = 0; isec < nsec; isec++ ) {
        Int_t secTrkId = fLastStackSize + isec; 
        // Tell secondary track about siblings
        PTSimParticle* secParticle = const_cast<PTSimParticle*>
                                   (fPTSimStack->GetParticle(secTrkId));
        for ( Int_t isec2 = 0; isec2 < nsec; isec2++ ) {
          if ( isec2 != isec ) {
            Int_t secTrkId2 = fLastStackSize + isec2;
            secParticle -> AddSibling(fPTSimStack->GetParticle(secTrkId2));
          }
        }
      }
    }
  }
    
  if ( fPTSimStack->GetStdHepSaveSibling() ) 
                            fLastStackSize = fPTSimStack->GetNtrack();

  // Flag to stop tracking if entering Mars from Liner
  if ( IsExiting() ) {
    MSG("PTSim",Msg::kDebug) 
      << "Track is entering MARS from LINR, call StopTrack " << endl;
    gMC -> StopTrack();  
    return;
  }
  if ( IsNeutrino() ) {
    MSG("PTSim",Msg::kDebug)
      << "Track is neutrino, call StopTrack " << endl;
    gMC -> StopTrack();
    return;
  }
  
  if ( !(gMC->TrackCharge()) ) return; // no hits recorded for uncharged

  bool isSensitive = false;
  Int_t currentMedId = gMC->CurrentMedium();
  for ( int imed = 0; imed < fNSensitiveMedId; imed++ ) {
    if ( fSensitiveMedId[imed] == currentMedId ) {
      isSensitive = true;
      MSG("PTSim",Msg::kVerbose) << "Sensitive volume" << endl;
      break;    
    }
  }
  
  // Only for sensitive volumes record hits
  if ( !isSensitive ) return;

  // Check to see if particle IsTrackStop and yet PTSimHit has not
  // been initialized.  This occurs in VMC w/G4 which double steps
  // some of the time when track has disappeared, e.g. to Decay or to
  // Positron annihilation.  The first call to Stepping will have
  // track status IsTrackDisappeared, at which point the PTSimHit
  // is closed.  The second call will have IsTrackStop.  The behavior
  // is inconsistent, and sometimes Stepping is only called once.
  // On the other hand, VMC w/G3 always only calls Stepping once,
  // with either status IsTrackDisappeared for decays, positron 
  // annhilation, etc. or IsTrackStop for energy below threshold. 
  if ( gMC -> IsTrackStop() && fMCHit.GetTrkIndex() < 0 ) return;
  
  if ( gMC -> IsTrackExiting() ) {
    // Don't trust gGeoManager->GetCurrentNode() directly for exiting particles
    // because TGeant4 doesn't set this properly before Stepping()
    fCurrentPath = gGeoManager->GetPath();
    gGeoManager->cd(gMC->CurrentVolPath());
  }

  // These should all be scint strip pstyrene insert nodes inside coex
  // shell, ask for parent GeoStripNode
  GeoStripNode* stp = dynamic_cast<GeoStripNode*>(gGeoManager->GetMother());
  if ( !stp ) {
    MSG("PTSim",Msg::kFatal) << "Stepping claims to be in sensitive volume "
                             << " but not strip node! abort." << endl;
    abort();
  }
  PlexStripEndId seid = stp -> GetSEId();

  // convert global position coordinates to coordinates local to strip
  TLorentzVector pos;
  gMC->TrackPosition(pos);
  TLorentzVector mom;
  gMC->TrackMomentum(mom);
  Double_t gxyz[3] = {pos.X(),pos.Y(),pos.Z()};
  Double_t lxyz[3] = {0};
  // Use GeoNode::GlobalToLocal and not TGeoManager::MasterToLocal so 
  // that coordinates will be stored in strip local coordinates, even 
  // when strip is split...
  // Modified S. Kasahara 2006/07/11. The following clause about
  // GeoNode::GlobalToLocal altering the state is no longer true.
  // Comment out the  "gGeoManager->cd(currentpath)" statement 
  // since this is cpu-intensive.
  // // but need to be careful because GeoNode::GlobalToLocal will
  // //  alter the state (cd to strip node path and not psty insert).
  // //  Save current global path and cd back after conversion is done.
  //    const char* currentpath = gGeoManager->GetPath();
  stp -> GlobalToLocal(gxyz,lxyz);
  //   gGeoManager -> cd(currentpath);

  TLorentzVector localpos;
  localpos.SetXYZT(lxyz[0],lxyz[1],lxyz[2],pos.T());

  Int_t trackid = fPTSimStack->GetCurrentTrackNumber();
  if ( gMC -> IsTrackEntering() ) {
    fMCHit.Clear();
    fMCHit.Init(seid,gMC->TrackPid(),trackid,pos,localpos,mom);
  }
  else {
    // accumulate summed path length & energy deposition
    if ( fMCHit.GetTrkIndex() != trackid ) {
      MSG("PTSim",Msg::kError) << "Current track no " << trackid
         << " not equal to hit trk index " 
         << fMCHit.GetTrkIndex() << ". abort." << endl;
      abort();
    }
    fMCHit.AddELoss(gMC->Edep());
    fMCHit.AddStep(gMC->TrackStep());
    fMCHit.SetCurrentState(pos,localpos,mom);
  }
 
  if ( gMC -> IsTrackExiting() || gMC -> IsTrackStop()
                               || gMC -> IsTrackDisappeared() ) {
    if ( fMCHit.GetTrkIndex() != trackid ) {
      MSG("PTSim",Msg::kError) << "Current track no " << trackid
         << " not equal to hit trk index " 
         << fMCHit.GetTrkIndex() << ". abort" << endl;
      abort();
    }
    PTSimParticle* track = const_cast<PTSimParticle*>
                             (fPTSimStack->GetParticle(trackid));
    if ( !(track->HasHitAboveThresh()) && 
         fMCHit.GetELoss() >= fPTSimStack->GetStdHepHitThr() ) {
      track -> SetHasHitAboveThresh(true);
    }
    
    fMCHit.CreateDigiScintHit(*fHitArray);
    fMCHit.Clear(); // clear for next track (overkill but to be safe)
  }

  if ( gMC -> IsTrackExiting() ) {
    // Don't trust gGeoManager->GetCurrentNode() directly for exiting particles
    // because TGeant4 doesn't set this properly before Stepping()
    gGeoManager->cd(fCurrentPath.c_str());
  }

}


//_____________________________________________________________________________
bool PTSimApplication::IsExiting() {    
  // Private method called on each step to determine if track is exiting. 
  // Determined as track has previously exited liner and is now entering
  // mars volume.  The algorithm is taken from R. Hatcher's gminos. 
  //

  if ( gMC -> IsTrackExiting() ) {
    if ( std::string(gMC -> CurrentVolName() ) == "LINR" ) {
      fExitLiner = true;
    }
  }
  else if ( gMC -> IsTrackEntering() && fExitLiner ) {
    if ( std::string(gMC -> CurrentVolName() ) == "MARS" ) {
      return true;
    }
  }

  return false;
  
}

//_____________________________________________________________________________
bool PTSimApplication::IsNeutrino() {    
  // Private method called on each step to determine if track particle
  // type is neutrino.   

  Int_t abspdgId = TMath::Abs(gMC -> TrackPid());

  if ( abspdgId == 12 ) return true;  // nu_e, nu_e_bar
  else if ( abspdgId == 14 ) return true;  // nu_mu, nu_mu_bar
  else if ( abspdgId == 16 ) return true;  // nu_tau, nu_tau_bar
  else if ( abspdgId == 18 ) return true;  // nu'_tau, nu'_tau_bar
  
  return false;
  
}

//_____________________________________________________________________________
void PTSimApplication::PostTrack() {    
  // User actions after finishing of each track

  MSG("PTSim",Msg::kVerbose) << "PTSimApplication::PostTrack" << endl;

}

//_____________________________________________________________________________
void PTSimApplication::FinishPrimary() {    
  // User actions after finishing of a primary track

  MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishPrimary " 
                       << fSnarl << "/" << fEvent << "/" << fPrimary << endl;

}

//_____________________________________________________________________________
void PTSimApplication::FinishEvent() {    
  // User actions after finishing of an event

  MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishEvent " 
                           << fSnarl << "/" << fEvent << endl;

  // Fill stdhep array with all primaries and selected secondaries
  // Adjust digihit track id's starting with first hit in this event 
  // to reflect stdheparray indices
  fPTSimStack -> FillStdHepArray(fStdHepArray,fSnarl,fHitArray,fNHitBegEvt);
  fPTSimStack -> Reset();

} 

//_____________________________________________________________________________
void PTSimApplication::InitSnarl(Int_t snarl, TClonesArray* stdheparray,
                                 Int_t nevent) {
  //
  // Purpose:  Private method used to fill fEvtStdHepMap with events
  //           and associated primaries based on information in input
  //           stdheparray.  At the end of the method, the stdheparray
  //           contents will be erased in preparation for the transport.
  //           At the end of transport the stdheparray will be filled
  //           with both primaries & secondaries. 
  //
  // Arguments: snarl: current snarl number
  //            stdheparray: filled with primaries and secondaries. Ownership
  //                         remains that of the caller.
  //            nevent: number of events in snarl (default 1)
  //  
  // Return:  none.
  //


  // Extract random number seeds to print.
  const MCApplication& mcapp = MCApplication::Instance();
  std::string evtstr = " event.";
  if ( nevent > 1 ) evtstr = " events.";
  
  MSG("PTSim",Msg::kInfo) << "***** PTSimApplication::InitSnarl: Snarl " 
                          << snarl << " with " << nevent << evtstr.c_str() 
                          << " Starting seed " 
                          << mcapp.GetRandom()->GetSeed() << " *****" << endl;

  if ( fRkVLevel > 0 && fRockdEdXMin < 0. ) InitRockdEdXMin();

  ResetEvtStdHepMap();
  fSnarl = snarl; // current snarl
  fNSnarls++;  // total number of snarls
  fEvent=-1;
  
  fStdHepArray = stdheparray; // reference to stdheparray, not owned.
  
  if ( nevent <= 0 ) return; // nothing to build

  Int_t nstdhep = fStdHepArray->GetEntriesFast();

  Int_t mcindex = -1;
  Int_t nprim = 0;
  Int_t idx0 = 0;
  bool  prevchild = true;
  std::map<int,int> indexmap;

  for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
    TParticle* simstdhep = dynamic_cast<TParticle*>(fStdHepArray->At(istd));

    if ( simstdhep->GetMother(0) == -1 && simstdhep->GetMother(1) == -1 ) {
      // primary, check to make sure its new event
      nprim++;
      if ( nprim > 2 || prevchild ) {
        if ( mcindex != -1 ) {

          // Clean up previous event stdhep indices
          std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
          Int_t nsize = particlelist.size();
          for ( int ip = 0; ip < nsize; ip++ ) {
            TParticle* particle = particlelist[ip];
            if ( particle->GetFirstMother() != -1 ) {
              // All mothers should be in index map
              std::map<int,int>::iterator inditr 
                                = indexmap.find(particle->GetFirstMother());
              if (inditr == indexmap.end() ) {
                MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
                                         << endl;
                abort();
              }
              particle->SetFirstMother(inditr->second);
            }
            if ( particle->GetSecondMother() != -1 ) {
              // All mothers should be in index map
              std::map<int,int>::iterator inditr 
                                = indexmap.find(particle->GetSecondMother());
              if (inditr == indexmap.end() ) {
                MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
                                         << endl;
                abort();
              }
              particle->SetLastMother(inditr->second);
            }
            if ( particle->GetFirstDaughter() != -1 ) {
              // check for first daughter in map
              std::map<int,int>::iterator inditr 
                                = indexmap.find(particle->GetFirstDaughter());
              if (inditr != indexmap.end() )
                particle->SetFirstDaughter(inditr->second);
              else
                particle->SetFirstDaughter(-1);
            }             
            if ( particle->GetLastDaughter() != -1 ) {
              // check for last daughter in map
              std::map<int,int>::iterator inditr 
                               = indexmap.find(particle->GetLastDaughter());
              if (inditr != indexmap.end() )
                particle->SetLastDaughter(inditr->second);
              else
                particle->SetLastDaughter(-1);
            }
          }
          indexmap.clear();
        }
      
        mcindex++;
        nprim = 1;
        idx0 = istd;
        if ( mcindex >= nevent ) {
          MSG("PTSim",Msg::kError)
             << "PTSimModule::BuildEvtStdHepMap:\nBreakdown in procedure "
             << "to match stdhep to associated mc entry.\n"
             << "Unable to process snarl.  Abort." << endl;
          abort();
        }
      }
    
      prevchild = false;
    }
    
    else {
      prevchild = true;
      nprim = 0;
    }
    

    Int_t statuscode = simstdhep->GetStatusCode();
    bool isprimary = false;
    if ( statuscode < 200 || statuscode == 999 ) {
      // Primary
      isprimary = true;
    }
    else if ( statuscode == 1205 ) {
      // May also be primary if status code of children is < 200
      TParticle* child = dynamic_cast<TParticle*>
                          (fStdHepArray->At(simstdhep->GetFirstDaughter()));
      if ( child->GetStatusCode() < 200 ) isprimary = true;
    }
    if ( isprimary ) {
      // Parent/child indices will be cleaned up at end of this event
      // Push primaries only
      std::vector<TParticle*>& stdheplist = fEvtStdHepMap[mcindex];
      TParticle* stdhepcopy = new TParticle(*simstdhep);
      // Convert to cm for use in simulation
      stdhepcopy -> SetProductionVertex(stdhepcopy->Vx()*100.,
                                        stdhepcopy->Vy()*100.,
                                        stdhepcopy->Vz()*100.,
                                        stdhepcopy->T());
      // Adjust pdg code if ion to be pdg2006 standard if not already done
      if ( UtilPDG::isIon(stdhepcopy -> GetPdgCode()) ) {
        int pdg2006 = UtilPDG::stdIonNumber(stdhepcopy->GetPdgCode(),
                                            UtilPDG::kPDG2006);
        stdhepcopy -> SetPdgCode(pdg2006);
      }
 
      indexmap.insert(std::make_pair(istd,stdheplist.size()));
      stdheplist.push_back(stdhepcopy);
    }
    
  }

  // At end, check grand total
  if ( mcindex != nevent - 1  ) {
    MSG("PTSim",Msg::kError)
      << "PTSimModule::InitSnarl:"
      << "\nBreakdown in procedure to match stdhep to "
      << "associated mc entry.\nFound " << mcindex + 1 << " primaries "
      << "in stdhep array, but nevent = " << nevent << ".\n"
      << "Unable to process snarl.  Abort." << endl;
    abort();
  }


  // Clean up last event stdhep indices
  std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
  Int_t nsize = particlelist.size();
  for ( int ip = 0; ip < nsize; ip++ ) {
    TParticle* particle = particlelist[ip];
    if ( particle->GetFirstMother() != -1 ) {
      // All mothers should be in index map
      std::map<int,int>::iterator inditr 
                 = indexmap.find(particle->GetFirstMother());
      if (inditr == indexmap.end() ) {
        MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
        abort();
      }
      particle->SetFirstMother(inditr->second);
    }
    if ( particle->GetSecondMother() != -1 ) {
      // All mothers should be in index map
      std::map<int,int>::iterator inditr 
                 = indexmap.find(particle->GetSecondMother());
      if (inditr == indexmap.end() ) {
        MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
        abort();
      }
      particle->SetLastMother(inditr->second);
    }
    if ( particle->GetFirstDaughter() != -1 ) {
      // check for daughter in map
      std::map<int,int>::iterator inditr 
                  = indexmap.find(particle->GetFirstDaughter());
      if (inditr != indexmap.end() )
        particle->SetFirstDaughter(inditr->second);
      else
        particle->SetFirstDaughter(-1);
    }             
    if ( particle->GetLastDaughter() != -1 ) {
      // check for daughter in map
      std::map<int,int>::iterator inditr 
                 = indexmap.find(particle->GetLastDaughter());
      if (inditr != indexmap.end() )
        particle->SetLastDaughter(inditr->second);
      else
        particle->SetLastDaughter(-1);
     }
   }
   
   indexmap.clear();

   // Erase stdheparray contents, pointers to primary TParticles
   // now contained in map
   fStdHepArray -> Clear("C");

   return;

}

//_____________________________________________________________________________
void PTSimApplication::ResetEvtStdHepMap() {

  EvtStdHepMap::reverse_iterator ritr;
  for ( ritr = fEvtStdHepMap.rbegin(); ritr != fEvtStdHepMap.rend(); ritr++ ) {
    std::vector<TParticle*>& particlelist = ritr->second;
    for ( unsigned int ip = 0; ip < particlelist.size(); ip++ ) {
      TParticle* particle = particlelist[ip];
      if ( particle ) delete particle; particle = 0;
    }
    particlelist.clear();
  }
  fEvtStdHepMap.clear();
  
}

//_____________________________________________________________________________
bool PTSimApplication::IsRockVetoed() {
  // Private method called on each primary to determine if track is in rock,
  // is charged, and fails to satisfy minimum energy threshold required to
  // reach detector cavern.  Return true if track is "Vetoed."

  bool isRockVetoed = false;
  
  // No veto'ing for uncharged particles
  if ( !(gMC->TrackCharge()) ) return isRockVetoed;

  Int_t copyNo = 0;
  if ( gMC -> CurrentVolID(copyNo) 
       == gGeoManager->GetMasterVolume()->GetNumber() ) {
    // In rock
    TGeoNode* node = gGeoManager->GetMasterVolume()->GetNode("LINR_1");
    TLorentzVector pos;
    gMC -> TrackPosition(pos);
    Double_t gxyz[3] = {pos.X(),pos.Y(),pos.Z()};
    Double_t safety = node -> Safety(gxyz,false);
    if ( safety > fRkVMinDistInCM ) {
      TLorentzVector mom;
      gMC->TrackMomentum(mom);
      if ( safety > fRkVFactor*mom.E()/fRockdEdXMin ) isRockVetoed = true;
    }
  }
  
  return isRockVetoed;
  
}

//_____________________________________________________________________________
void PTSimApplication::InitRockdEdXMin() {
  // Private method to initialize fRockdEdXMin with the minimum dE/dX from 
  // LOSS tables for MARS rock material for mu-. Units of fRockdEdXMin are 
  // GeV/cm.  This requires the use of TGeant3 specific code to do it right, 
  // which is to check the LOSS tables.  For non-TGeant3 concrete MC, will use 
  // a nominal dE/dX value of 1.2e-3 GeV/(g/cm^2) * density of the rock. 


  Double_t nominaldEdXMin = 1.2e-3;  // GeV/(g/cm^2)

  TGeant3* geant3 = dynamic_cast<TGeant3*>(gMC);
  if ( !geant3 ) {
    MSG("PTSim",Msg::kWarning) 
      << "No support for non-TGeant3 concrete VMC in InitRockdEdXMin\n"
      << "to determine minimum dE/dX from LOSS tables.  Will use nominal "
      << "value of " << nominaldEdXMin <<  " GeV/(g/cm^2) * the density "
      << " of the rock." << endl;
    fRockdEdXMin = nominaldEdXMin 
             * (gGeoManager->GetVolume("MARS")->GetMaterial()->GetDensity());
  }
  else {
    // TGeant3 allows access to LOSS tables, so do it the right way
    
    // Determine kinetic energy range of interest using gcmulo common block to 
    // extract range for which LOSS table has been calculated.  Nominal range 
    // is 10^-5 GeV to 10^4 GeV.  
    Float_t log10KEmin = log10(geant3 -> Gcmulo() -> ekmin);
    log10KEmin = TMath::Max(log10KEmin,(Float_t)-5.); // don't go below 10 keV
    Float_t log10KEmax = log10(geant3 -> Gcmulo() -> ekmax);

    // Now prepare to use Gftmat to extract dE/dX

    // Set up 1000 input KE bins (tkin) at which dE/dX will be extracted from 
    // table (value).
    const Int_t nlog10KEbin = 1000;
    Float_t* tkin = new Float_t[nlog10KEbin];
    Float_t* value = new Float_t[nlog10KEbin];

    // Initialize output array returned by Gftmat, contains five energy 
    // thresholds for material
    Float_t pcut[5] = {0};

    // Material index and particle Id, const_cast because 
    // TGeoMaterial::GetIndex() isn't const 
    TString volName = "MARS";
    TString matName;  // material Name
    Int_t imate = 0;  // material Index
    Double_t a, z, dens, radl, inter;
    TArrayD par;
    if ( !(geant3 
           -> GetMaterial(volName,matName,imate,a,z,dens,radl,inter,par)) ) {
      MSG("PTSim",Msg::kError) << "Failed to retrieve material for volume " 
                               << volName.Data() << ". Abort." << endl;
      abort();
    }
    MSG("PTSim",Msg::kDebug) 
      << "InitRockdEdXMin: GetMaterial for volume MARS retrieved:\n"
      << " imate " << imate << " a " << a << " z " << z << " dens " << dens 
      << " name " << matName.Data() << endl;
      
    Int_t ipart = 6;  // mu-
    
    char mechname[] = "LOSS";

    Float_t de = (log10KEmax - log10KEmin)/(nlog10KEbin);
    for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
      Float_t log10tkin = log10KEmin + Float_t(ie)*de;
      tkin[ie] = pow((Float_t)10., log10tkin);
      value[ie] = 0;
    }
    // Gftmat input:
    //    imate - material number
    //    ipart - particle number
    //    chmeca - name of mechanism bank to be retrieved
    //    kdim - dimension of the arrays tkin and value
    //    tkin - array of kinetic energies in GeV
    // output:
    //    value - array of energy losses in MeV/cm or macroscopic cross 
    //            sections in 1/cm
    //    pcut - array containing the 5 energy thresholds for the materials
    //           in GeV (CUTGAM, CUTELE, CUTNEU, CUTHAD, CUTMUO)
    //    ixst - return code ( = 1 if filled, = 0 otherwise)
    //  
    Int_t ixst = 0;
    geant3 -> Gftmat(imate,ipart,mechname,nlog10KEbin,tkin,value,pcut,ixst);
    
    if ( ixst != 1 ) {
      MSG("PTSim",Msg::kError) << "Gftmat call to retrieve dE/dX for matId " 
                               << imate << " over log10(KE) range "
                               << log10KEmin << " to " << log10KEmax 
                               << "\nfailed when initializing energy "
                               << "threshold cuts.  Abort." << endl;
      abort();
    }    

    // Now loop over value, looking for min value.  This is min dE/dX
    Float_t minValue = 1.e10;
    for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
      minValue = TMath::Min(value[ie],minValue);
    }

    fRockdEdXMin = minValue*1.e-3; // convert to GeV/cm

    delete [] tkin;
    delete [] value;
    
  }
  
  MSG("PTSim",Msg::kInfo) 
    << "Initialize rock minimum dE/dX for application of veto cuts to: "
    << fRockdEdXMin << " (GeV/cm)" << endl;

  return;
  
}

  

---