Stack.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#include "DetectorsBase/Stack.h"
#include "DetectorsBase/Detector.h"
#include <DetectorsBase/VMCSeederService.h>
#include "DetectorsCommonDataFormats/DetID.h"
#include "SimulationDataFormat/MCTrack.h"
 
#include "FairDetector.h"      // for FairDetector
#include <fairlogger/Logger.h> // for FairLogger
#include "FairRootManager.h"
#include "SimulationDataFormat/BaseHits.h"
#include "SimulationDataFormat/StackParam.h"
#include "CommonUtils/ConfigurationMacroHelper.h"
 
#include "TLorentzVector.h" // for TLorentzVector
#include "TParticle.h"      // for TParticle
#include "TRefArray.h"      // for TRefArray
#include "TVirtualMC.h"     // for VMC
#include "TMCProcess.h"     // for VMC Particle Production Process
 
#include <algorithm>
#include <cassert>
#include <cstddef> // for NULL
#include <cmath>
#include <unordered_map>
#include "DetectorsBase/SimRandom.h"
 
using std::cout;
using std::endl;
using std::pair;
using namespace o2::data;
 
namespace
{
// a TParticle copied from a per-thread prototype per PDG code, avoiding the ROOT-locked PDG lookup
TParticle makeParticle(int pdg, int status, int mother1, int mother2, int daughter1, int daughter2,
                       double px, double py, double pz, double e, double vx, double vy, double vz, double t)
{
  thread_local std::unordered_map<int, TParticle> prototypes;
  auto it = prototypes.find(pdg);
  if (it == prototypes.end()) {
    TParticle p(pdg, status, mother1, mother2, daughter1, daughter2, px, py, pz, e, vx, vy, vz, t);
    if (p.GetPDG(1)) { // unknown codes may still be added to TDatabasePDG later
      prototypes.emplace(pdg, p);
    }
    return p;
  }
  TParticle p(it->second); // the copy keeps the PDG record and the mass
  p.SetStatusCode(status);
  p.SetFirstMother(mother1);
  p.SetLastMother(mother2);
  p.SetFirstDaughter(daughter1);
  p.SetLastDaughter(daughter2);
  p.SetMomentum(px, py, pz, e);
  p.SetProductionVertex(vx, vy, vz, t);
  return p;
}
} // namespace
 
// small helper function to append to vector at arbitrary position
template <typename T, typename I>
void insertInVector(std::vector<T>& v, I index, T e)
{
  auto currentsize = v.size();
  if (index >= currentsize) {
    const auto newsize = std::max(index + 1, (I)(1 + currentsize * 1.2));
    v.resize(newsize, T(-1));
  }
  // new size must at least be as large as index
  // assert(index < v.size());
  v[index] = e;
}
 
Stack::Stack(Int_t size)
  : FairGenericStack(),
    mStack(),
    mParticles(),
    mTracks(new std::vector<o2::MCTrack>),
    mTrackIDtoParticlesEntry(1000000, -1),
    mIndexMap(),
    mIndexOfCurrentTrack(-1),
    mIndexOfCurrentPrimary(-1),
    mNumberOfPrimaryParticles(0),
    mNumberOfEntriesInParticles(0),
    mNumberOfEntriesInTracks(0),
    mNumberOfPrimariesforTracking(0),
    mNumberOfPrimariesPopped(0),
    mStoreMothers(kTRUE),
    mStoreSecondaries(kTRUE),
    mMinHits(1),
    mEnergyCut(0.),
    mTrackRefs(new std::vector<o2::TrackReference>),
    mIsG4Like(false)
{
  auto vmc = TVirtualMC::GetMC();
  if (vmc) {
    mIsG4Like = !(vmc->SecondariesAreOrdered());
  }
 
  auto& param = o2::sim::StackParam::Instance();
  TransportFcn transportPrimary;
  if (param.transportPrimary.compare("none") == 0) {
    transportPrimary = [](const TParticle& p, const std::vector<TParticle>& particles) {
      return false;
    };
  } else if (param.transportPrimary.compare("all") == 0) {
    transportPrimary = [](const TParticle& p, const std::vector<TParticle>& particles) {
      return true;
    };
  } else if (param.transportPrimary.compare("barrel") == 0) {
    transportPrimary = [](const TParticle& p, const std::vector<TParticle>& particles) {
      return (std::fabs(p.Eta()) < 2.0);
    };
  } else if (param.transportPrimary.compare("external") == 0) {
    transportPrimary = o2::conf::GetFromMacro<o2::data::Stack::TransportFcn>(param.transportPrimaryFileName,
                                                                             param.transportPrimaryFuncName,
                                                                             "o2::data::Stack::TransportFcn", "stack_transport_primary");
    if (!transportPrimary) {
      LOG(fatal) << "Failed to retrieve external \'transportPrimary\' function: problem with configuration ";
    }
    LOG(info) << "Successfully retrieve external \'transportPrimary\' frunction: " << param.transportPrimaryFileName;
  } else {
    LOG(fatal) << "unsupported \'trasportPrimary\' mode: " << param.transportPrimary;
  }
 
  if (param.transportPrimaryInvert) {
    mTransportPrimary = [transportPrimary](const TParticle& p, const std::vector<TParticle>& particles) { return !transportPrimary(p, particles); };
  } else {
    mTransportPrimary = transportPrimary;
  }
}
 
Stack::Stack(const Stack& rhs)
  : FairGenericStack(rhs),
    mStack(),
    mParticles(),
    mTrackIDtoParticlesEntry(rhs.mTrackIDtoParticlesEntry.size(), -1),
    mTracks(nullptr),
    mIndexMap(),
    mIndexOfCurrentTrack(-1),
    mIndexOfCurrentPrimary(-1),
    mNumberOfPrimaryParticles(0),
    mNumberOfEntriesInParticles(0),
    mNumberOfEntriesInTracks(0),
    mNumberOfPrimariesforTracking(0),
    mNumberOfPrimariesPopped(0),
    mStoreMothers(rhs.mStoreMothers),
    mStoreSecondaries(rhs.mStoreSecondaries),
    mPruneKinematics(rhs.mPruneKinematics),
    mMinHits(rhs.mMinHits),
    mEnergyCut(rhs.mEnergyCut),
    mCleanupThreshold(rhs.mCleanupThreshold),
    mIsG4Like(rhs.mIsG4Like),
    mIsExternalMode(rhs.mIsExternalMode),
    mDoTrackSeeding(rhs.mDoTrackSeeding),
    mTransportPrimary(rhs.mTransportPrimary),
    mTrackRefs(new std::vector<o2::TrackReference>)
{
  LOG(debug) << "copy constructor called";
  mTracks = new std::vector<MCTrack>();
}
 
Stack::~Stack()
{
  if (mTracks) {
    delete mTracks;
  }
}
 
Stack& Stack::operator=(const Stack& rhs)
{
  LOG(fatal) << "operator= called";
  // check assignment to self
  if (this == &rhs) {
    return *this;
  }
 
  // base class assignment
  FairGenericStack::operator=(rhs);
  mTracks = new std::vector<MCTrack>(rhs.mTracks->size());
  mIndexOfCurrentTrack = -1;
  mIndexOfCurrentPrimary = -1;
  mNumberOfPrimaryParticles = 0;
  mNumberOfEntriesInParticles = 0;
  mNumberOfEntriesInTracks = 0;
  mNumberOfPrimariesforTracking = 0;
  mNumberOfPrimariesPopped = 0;
  mStoreMothers = rhs.mStoreMothers;
  mStoreSecondaries = rhs.mStoreSecondaries;
  mMinHits = rhs.mMinHits;
  mEnergyCut = rhs.mEnergyCut;
  mIsG4Like = rhs.mIsG4Like;
 
  return *this;
}
 
void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px, Double_t py, Double_t pz, Double_t e,
                      Double_t vx, Double_t vy, Double_t vz, Double_t time, Double_t polx, Double_t poly, Double_t polz,
                      TMCProcess proc, Int_t& ntr, Double_t weight, Int_t is)
{
  PushTrack(toBeDone, parentId, pdgCode, px, py, pz, e, vx, vy, vz, time, polx, poly, polz, proc, ntr, weight, is, -1);
}
 
void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px, Double_t py, Double_t pz, Double_t e,
                      Double_t vx, Double_t vy, Double_t vz, Double_t time, Double_t polx, Double_t poly, Double_t polz,
                      TMCProcess proc, Int_t& ntr, Double_t weight, Int_t is, Int_t secondparentId)
{
  PushTrack(toBeDone, parentId, pdgCode, px, py, pz, e, vx, vy, vz, time, polx, poly, polz, proc, ntr, weight, is, secondparentId, -1, -1, proc);
}
 
void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px, Double_t py, Double_t pz, Double_t e,
                      Double_t vx, Double_t vy, Double_t vz, Double_t time, Double_t polx, Double_t poly, Double_t polz,
                      TMCProcess proc, Int_t& ntr, Double_t weight, Int_t is, Int_t secondparentId, Int_t daughter1Id, Int_t daughter2Id,
                      TMCProcess proc2)
{
  //  printf("Pushing %s toBeDone %5d parentId %5d pdgCode %5d is %5d entries %5d \n",
  //   proc == kPPrimary ? "Primary:   " : "Secondary: ",
  //   toBeDone, parentId, pdgCode, is, mNumberOfEntriesInParticles);
 
  //
  // This method is called
  //
  // - during serial   simulation directly from the event generator
  // - during parallel simulation to fill the stack of the primary generator device
  // - in all cases to push a secondary particle
  //
  //
  // Create new TParticle and add it to the TParticle array
 
  Int_t trackId = mNumberOfEntriesInParticles;
  // Set track variable
  ntr = trackId;
  //  Int_t daughter1Id = -1;
  //  Int_t daughter2Id = -1;
  Int_t iStatus = (proc == kPPrimary) ? is : trackId;
  TParticle p = makeParticle(pdgCode, iStatus, parentId, secondparentId, daughter1Id, daughter2Id, px, py, pz, e, vx, vy, vz, time);
  p.SetPolarisation(polx, poly, polz);
  p.SetWeight(weight);
  p.SetUniqueID(proc);                                           // using the unique ID to transfer process ID
  p.SetBit(ParticleStatus::kPrimary, proc == kPPrimary ? 1 : 0); // set primary bit
  p.SetBit(ParticleStatus::kToBeDone, toBeDone == 1 ? 1 : 0);    // set to be done bit
  mNumberOfEntriesInParticles++;
 
  insertInVector(mTrackIDtoParticlesEntry, trackId, (int)(mParticles.size()));
 
  handleTransportPrimary(p); // handle selective transport of primary particles
 
  // Push particle on the stack if toBeDone is set
  if (p.TestBit(ParticleStatus::kPrimary)) {
    // This is a particle from the primary particle generator
    //
    // SetBit is used to pass information about the primary particle to the stack during transport.
    // Sime particles have already decayed or are partons from a shower. They are needed for the
    // event history in the stack, but not for transport.
    //
 
    // primary particles might have been pushed with a second creation process
    // in case we pushed a secondary track of a previous simulation to be continued.
    // We save therefore in the UniqueID the correct process
    // while the particle will still be treated as a primary given its bit settings
    p.SetUniqueID(proc2);
 
    mIndexMap[trackId] = trackId;
    p.SetBit(ParticleStatus::kKeep, 1);
    if (p.TestBit(ParticleStatus::kToBeDone)) {
      mNumberOfPrimariesforTracking++;
    }
    mNumberOfPrimaryParticles++;
    mPrimaryParticles.push_back(p);
    mTracks->emplace_back(p);
  } else {
    mParticles.emplace_back(p);
    mCurrentParticle0 = p;
  }
  mStack.push(p);
}
 
void Stack::handleTransportPrimary(TParticle& p)
{
  // this function tests whether we really want to transport
  // this particle and sets the relevant bits accordingly
 
  if (!p.TestBit(ParticleStatus::kToBeDone) || !p.TestBit(ParticleStatus::kPrimary)) {
    return;
  }
 
  if (!mTransportPrimary(p, mPrimaryParticles)) {
    p.SetBit(ParticleStatus::kToBeDone, 0);
    p.SetBit(ParticleStatus::kInhibited, 1);
  }
}
 
void Stack::PushTrack(int toBeDone, TParticle& p)
{
  //  printf("stack -> Pushing Primary toBeDone %5d %5d parentId %5d pdgCode %5d is %5d entries %5d \n", toBeDone, p.TestBit(ParticleStatus::kToBeDone), p.GetFirstMother(), p.GetPdgCode(), p.GetStatusCode(),  mNumberOfEntriesInParticles);
 
  // This method is called
  //
  // - during parallel simulation to push primary particles (called by the stack itself)
  if (p.TestBit(ParticleStatus::kPrimary)) {
    // one to one mapping for primaries
    mIndexMap[mNumberOfPrimaryParticles] = mNumberOfPrimaryParticles;
    mNumberOfPrimaryParticles++;
    mPrimaryParticles.push_back(p);
    // Push particle on the stack
    if (p.TestBit(ParticleStatus::kToBeDone)) {
      mNumberOfPrimariesforTracking++;
    }
    mStack.push(p);
    mTracks->emplace_back(p);
  }
}
 
void Stack::SetCurrentTrack(Int_t iTrack)
{
  mIndexOfCurrentTrack = iTrack;
  if (iTrack < mPrimaryParticles.size()) {
    auto& p = mPrimaryParticles[iTrack];
    mCurrentParticle = p;
    mIndexOfCurrentPrimary = iTrack;
  } else {
    mCurrentParticle = mCurrentParticle0;
  }
}
 
// calculates a hash based on particle properties
// hash may serve as seed for this track
ULong_t getHash(TParticle const& p)
{
  auto asLong = [](double x) {
    return (ULong_t) * (reinterpret_cast<ULong_t*>(&x));
  };
 
  ULong_t hash;
  o2::MCTrackT<double> track(p);
 
  hash = asLong(track.GetStartVertexCoordinatesX());
  hash ^= asLong(track.GetStartVertexCoordinatesY());
  hash ^= asLong(track.GetStartVertexCoordinatesZ());
  hash ^= asLong(track.GetStartVertexCoordinatesT());
  hash ^= asLong(track.GetStartVertexMomentumX());
  hash ^= asLong(track.GetStartVertexMomentumY());
  hash ^= asLong(track.GetStartVertexMomentumZ());
  hash += (ULong_t)track.GetPdgCode();
  return hash;
}
 
TParticle* Stack::PopNextTrack(Int_t& iTrack)
{
 
  // This functions is mainly used by Geant3?
  Int_t nprod = (int)(mParticles.size());
 
  // If end of stack: Return empty pointer
  if (mStack.empty()) {
    iTrack = -1;
    return nullptr;
  }
  Bool_t found = kFALSE;
 
  TParticle* nextParticle = nullptr;
  while (!found && !mStack.empty()) {
    // get next particle from stack
    mCurrentParticle = mStack.top();
    // remove particle from the top
    mStack.pop();
    // test if primary to be transported
    if (mCurrentParticle.TestBit(ParticleStatus::kToBeDone)) {
      if (mCurrentParticle.TestBit(ParticleStatus::kPrimary)) {
        // particle is primary and needs to be tracked -> indicates that previous particle finished
        mNumberOfPrimariesPopped++;
        mIndexOfCurrentPrimary = mStack.size();
        mIndexOfCurrentTrack = mIndexOfCurrentPrimary;
      } else {
        mIndexOfCurrentTrack = mCurrentParticle.GetStatusCode();
      }
      iTrack = mIndexOfCurrentTrack;
      if (mDoTrackSeeding) {
        auto hash = getHash(mCurrentParticle);
        // LOG(info) << "SEEDING NEW TRACK USING HASH" << hash;
        // init seed per track
        o2::base::simRandom()->SetSeed(hash);
        o2::base::VMCSeederService::instance().setSeed();
        // NOTE: THE BETTER PLACE WOULD BE IN PRETRACK HOOK BUT THIS DOES NOT SEEM TO WORK
        // WORKS ONLY WITH G3 SINCE G4 DOES NOT CALL THIS FUNCTION
      } // .trackSeed ?
      nextParticle = &mCurrentParticle;
      found = kTRUE;
    } else {
      iTrack = -1;
      nextParticle = nullptr;
    }
  } // while
  return nextParticle;
}
 
TParticle* Stack::PopPrimaryForTracking(Int_t iPrim)
{
  // This function is used by Geant4 to setup their own internal stack
  //  printf("PopPrimary for tracking %5d %5d \n", iPrim, mNumberOfPrimaryParticles);
  // Remark: Contrary to what the interface name is suggesting
  // this is not a pop operation (but rather a get)
 
  // Test for index
  if (iPrim < 0 || iPrim >= mNumberOfPrimaryParticles) {
    LOG(fatal) << "Stack::PopPrimaryForTracking: Stack: Primary index out of range! " << iPrim << " ";
    return nullptr;
  }
  // Return the iPrim-th TParticle from the fParticle array. This should be
  // a primary.
  auto particle = &mPrimaryParticles[iPrim];
  if (particle->TestBit(ParticleStatus::kToBeDone)) {
    return particle;
  } else {
    return nullptr;
  }
}
 
void Stack::updateEventStats()
{
  if (mMCEventStats) {
    mMCEventStats->setNHits(mHitCounter);
    mMCEventStats->setNTransportedTracks(mNumberOfEntriesInParticles);
    mMCEventStats->setNKeptTracks(mTracks->size());
  }
}
 
void Stack::FillTrackArray()
{
  LOG(info) << "Stack: " << mTracks->size() << " out of " << mNumberOfEntriesInParticles << " stored \n";
}
 
void Stack::FinishPrimary()
{
  // Here transport of a primary and all its secondaries is finished
  // we can do some cleanup of the memory structures
  mPrimariesDone++;
  LOG(debug) << "Finish primary hook " << mPrimariesDone;
  // preserve particles and theire ancestors that produced hits
 
  auto selected = selectTracks();
 
  // loop over current particle buffer
  // - build index map indicesKept
  // - update mother index information
  int indexOld = 0;
  int indexNew = 0;
  int indexoffset = mTracks->size();
  int neglected = 0;
  std::vector<int> indicesKept((int)(mParticles.size()));
  std::vector<MCTrack> tmpTracks;
  Int_t ic = 0;
 
  // mTrackIDtoParticlesEntry
  // trackID to mTrack -> index in mParticles
  //
  // indicesKept
  // index in mParticles -> index in mTrack - highWaterMark
  //
  // mReorderIndices
  // old (mTrack-highWaterMark) -> new (mTrack-highWaterMark)
  //
 
  for (auto& particle : mParticles) {
    if (particle.getStore() || !mPruneKinematics) {
      // map the global track index to the new persistent index
      auto imother = particle.getMotherTrackId();
      // here mother is relative to the mParticles array
      if (imother >= 0) {
        // daughter of a secondary: obtain index from lookup table
        imother = indicesKept[imother];
        particle.SetMotherTrackId(imother);
      }
      // at this point we have the correct mother index in mParticles or
      // a negative one which is a pointer to a primary
      tmpTracks.emplace_back(particle);
      indicesKept[indexOld] = indexNew;
      indexNew++;
    } else {
      indicesKept[indexOld] = -1;
      neglected++;
    }
    indexOld++;
    mTracksDone++;
  }
  Int_t ntr = (int)(tmpTracks.size());
  std::vector<int> reOrderedIndices(ntr);
  std::vector<int> invreOrderedIndices(ntr);
  for (Int_t i = 0; i < ntr; i++) {
    invreOrderedIndices[i] = i;
    reOrderedIndices[i] = i;
  }
 
  if (mIsG4Like) {
    ReorderKine(tmpTracks, reOrderedIndices);
    for (Int_t i = 0; i < ntr; i++) {
      Int_t index = reOrderedIndices[i];
      invreOrderedIndices[index] = i;
    }
  }
  for (Int_t i = 0; i < ntr; i++) {
    Int_t index = reOrderedIndices[i];
    auto& particle = tmpTracks[index];
    Int_t imo = particle.getMotherTrackId();
    Int_t imo0 = imo;
    if (imo >= 0) {
      imo = invreOrderedIndices[imo];
    }
    imo += indexoffset;
    particle.SetMotherTrackId(imo);
    mTracks->emplace_back(particle);
    auto& mother = mTracks->at(imo);
    if (mother.getFirstDaughterTrackId() == -1) {
      mother.SetFirstDaughterTrackId((int)(mTracks->size()) - 1);
    }
    mother.SetLastDaughterTrackId((int)(mTracks->size()) - 1);
  }
 
  //
  // Update index map
  //
  Int_t imax = mNumberOfEntriesInParticles;
  Int_t imin = imax - mParticles.size();
  for (Int_t idTrack = imin; idTrack < imax; idTrack++) {
    Int_t index1 = mTrackIDtoParticlesEntry[idTrack];
    Int_t index2 = indicesKept[index1];
    if (index2 == -1) {
      continue;
    }
    Int_t index3 = (mIsG4Like) ? invreOrderedIndices[index2] : index2;
    mIndexMap[idTrack] = index3 + indexoffset;
  }
 
  // we can now clear the particles buffer!
  reOrderedIndices.clear();
  invreOrderedIndices.clear();
  mParticles.clear();
  tmpTracks.clear();
  mTransportedIDs.clear();
  mTrackIDtoParticlesEntry.clear();
  mIndexOfPrimaries.clear();
}
 
void Stack::UpdateTrackIndex(TRefArray* detList)
{
  // we can avoid any updating in case no tracks have been filtered out
  // check this like this
  if (mIndexMap.size() == 0) {
    LOG(info) << "No TrackIndex update necessary\n";
    return;
  }
 
  // we are getting the detectorlist from FairRoot as TRefArray
  // at each call, but this list never changes so we cache it here
  // as the right type to avoid repeated dynamic casts
  if (mActiveDetectors.size() == 0) {
    if (detList == nullptr) {
      LOG(fatal) << "No detList passed to Stack";
    }
    auto iter = detList->MakeIterator();
    while (auto det = iter->Next()) {
      auto o2det = dynamic_cast<o2::base::Detector*>(det);
      if (o2det) {
        mActiveDetectors.emplace_back(o2det);
      } else {
        LOG(info) << "Found nonconforming detector " << det->GetName();
      }
    }
  }
 
  LOG(debug) << "Stack::UpdateTrackIndex: Stack: Updating track indices...";
  Int_t nColl = 0;
 
  // update track references
  // use some caching since repeated trackIDs
  for (auto& ref : *mTrackRefs) {
    const auto id = ref.getTrackID();
    auto iter = mIndexMap.find(id);
    if (iter == mIndexMap.end()) {
      LOG(info) << "Invalid trackref ... needs to be rmoved \n";
      ref.setTrackID(-1);
    } else {
      ref.setTrackID(iter->second);
    }
  }
 
  // sort trackrefs according to new track index
  // then according to track length
  std::sort(mTrackRefs->begin(), mTrackRefs->end(), [](const o2::TrackReference& a, const o2::TrackReference& b) {
    if (a.getTrackID() == b.getTrackID()) {
      return a.getLength() < b.getLength();
    }
    return a.getTrackID() < b.getTrackID();
  });
 
  for (auto det : mActiveDetectors) {
    // update the track indices by delegating to specialized detector functions
    det->updateHitTrackIndices(mIndexMap);
  } // List of active detectors
 
  LOG(debug) << "Stack::UpdateTrackIndex: ...stack and " << nColl << " collections updated.";
}
 
void Stack::Reset()
{
 
  mIndexOfCurrentTrack = -1;
  mNumberOfPrimaryParticles = mNumberOfEntriesInParticles = mNumberOfEntriesInTracks = 0;
  while (!mStack.empty()) {
    mStack.pop();
  }
  mParticles.clear();
  mTracks->clear();
  if (!mIsExternalMode && (mPrimariesDone != mNumberOfPrimariesforTracking)) {
    LOG(fatal) << "Inconsistency in primary particles treated " << mPrimariesDone << " vs expected "
               << mNumberOfPrimariesforTracking << "\n(This points to a flaw in the stack logic)";
  }
  mPrimariesDone = 0;
  mNumberOfPrimariesforTracking = 0;
  mNumberOfPrimariesPopped = 0;
  mPrimaryParticles.clear();
  mTrackRefs->clear();
  mTrackIDtoParticlesEntry.clear();
  mIndexMap.clear();
  mHitCounter = 0;
}
 
void Stack::Register()
{
  FairRootManager::Instance()->RegisterAny("MCTrack", mTracks, kTRUE);
  FairRootManager::Instance()->RegisterAny("TrackRefs", mTrackRefs, kTRUE);
}
 
void Stack::Print(Int_t iVerbose) const
{
  cout << "-I- Stack: Number of primaries        = " << mNumberOfPrimaryParticles << endl;
  cout << "              Total number of particles  = " << mNumberOfEntriesInParticles << endl;
  cout << "              Number of tracks in output = " << mNumberOfEntriesInTracks << endl;
  if (iVerbose) {
    for (auto& track : *mTracks) {
      track.Print();
    }
  }
}
 
void Stack::Print(Option_t* option) const
{
  Int_t verbose = 0;
  if (option) {
    verbose = 1;
  }
  Print(verbose);
}
 
void Stack::addHit(int iDet)
{
  // translate detector ID to bit id
  const auto bitID = o2::base::Detector::getDetId2HitBitIndex()[iDet];
 
  if (mIndexOfCurrentTrack < mNumberOfPrimaryParticles) {
    auto& part = mTracks->at(mIndexOfCurrentTrack);
    part.setHit(bitID);
 
  } else {
    Int_t iTrack = mTrackIDtoParticlesEntry[mIndexOfCurrentTrack];
    auto& part = mParticles[iTrack];
    part.setHit(bitID);
  }
  mCurrentParticle.SetBit(ParticleStatus::kHasHits, 1);
  mHitCounter++;
}
void Stack::modifiesHit(int iDet)
{
  // flag only; hit mask, kHasHits and hit counter stay untouched
  if (mIndexOfCurrentTrack < mNumberOfPrimaryParticles) {
    mTracks->at(mIndexOfCurrentTrack).setContributesToAHit();
  } else {
    Int_t iTrack = mTrackIDtoParticlesEntry[mIndexOfCurrentTrack];
    mParticles[iTrack].setContributesToAHit();
  }
}
void Stack::addHit(int iDet, Int_t iTrack)
{
  mHitCounter++;
  auto& part = mParticles[iTrack];
  // fetch the bit encoding for hits
  const auto bitID = o2::base::Detector::getDetId2HitBitIndex()[iDet];
  part.setHit(bitID);
  mCurrentParticle.SetBit(ParticleStatus::kHasHits, 1);
}
 
Int_t Stack::GetCurrentParentTrackNumber() const
{
  TParticle* currentPart = GetCurrentTrack();
  if (currentPart) {
    return currentPart->GetFirstMother();
  } else {
    return -1;
  }
}
 
bool Stack::selectTracks()
{
  bool tracksdiscarded = false;
  // Check particles in the fParticle array
  int prim = -1; // counter how many primaries seen (mainly to constrain search in motherindex remapping)
  LOG(debug) << "Stack: Entering track selection on " << mParticles.size() << " tracks";
  for (auto& thisPart : mParticles) {
    Bool_t store = kTRUE;
    // Get track parameters
    Bool_t isPrimary = (thisPart.getProcess() == 0);
    Int_t iMother = thisPart.getMotherTrackId();
    auto& mother = mParticles[iMother];
    Bool_t motherIsPrimary = (iMother < mNumberOfPrimaryParticles);
 
    if (isPrimary) {
      prim++;
      // for primaries we are done quickly
      // in fact the primary has already been stored
      // so this should not happen
      store = kTRUE;
    } else {
      // for other particles we potentially need to correct the mother indices
      if (!motherIsPrimary) {
        // the mapping is from the index in the stack to the index in mParticles
        thisPart.SetMotherTrackId(mTrackIDtoParticlesEntry[iMother]);
      } else {
        // for a secondary which is a direct decendant of a primary use a negative index which will be restored later
        thisPart.SetMotherTrackId(iMother - (int)(mTracks->size()));
      }
      // no secondaries; also done
      if (!mStoreSecondaries) {
        store = kFALSE;
        tracksdiscarded = true;
      } else {
        // Calculate number of hits created by this track
        // Note: we only distinguish between no hit and more than 0 hits
        int nHits = thisPart.hasHits();
        // Check for cuts (store primaries in any case)
        if (nHits < mMinHits) {
          store = kFALSE;
          tracksdiscarded = true;
        }
        // only if we have non-trival energy cut
        if (mEnergyCut > 0.) {
          Double_t energy = thisPart.GetEnergy();
          Double_t mass = thisPart.GetMass();
          Double_t eKin = energy - mass;
          if (eKin < mEnergyCut) {
            store = kFALSE;
            tracksdiscarded = true;
          }
        }
        if (keepPhysics(thisPart)) {
          store = kTRUE;
          tracksdiscarded = false;
        }
      }
    }
 
    store = store || thisPart.getStore();
    thisPart.setStore(store);
  }
 
  // If flag is set, flag recursively mothers of selected tracks
  //
  if (mStoreMothers) {
    for (auto& particle : mParticles) {
      if (particle.getStore()) {
        Int_t iMother = particle.getMotherTrackId();
        while (iMother >= 0) {
          auto& mother = mParticles[iMother];
          mother.setStore(true);
          iMother = mother.getMotherTrackId();
        } // while mother
      }   // store ?
    }     // particle loop
  }
 
  return !tracksdiscarded;
}
 
bool Stack::isPrimary(const MCTrack& part)
{
  if (part.getProcess() == kPPrimary || part.getMotherTrackId() < 0) {
    return true;
  }
  return false;
}
 
bool Stack::isFromPrimaryDecayChain(const MCTrack& part)
{
  if (part.getProcess() != kPDecay) {
    return false;
  }
  auto imother = part.getMotherTrackId();
  auto mother = mParticles[imother];
  if (isPrimary(mother)) {
    return true;
  }
  return isFromPrimaryDecayChain(mother);
}
 
bool Stack::isFromPrimaryPairProduction(const MCTrack& part)
{
  if (part.getProcess() != kPPair) {
    return false;
  }
  auto imother = part.getMotherTrackId();
  auto mother = mParticles[imother];
  if (isPrimary(mother)) {
    return true;
  }
  return isFromPrimaryDecayChain(mother);
}
 
bool Stack::keepPhysics(const MCTrack& part)
{
  //
  // Some particles have to kept on the stack for reasons motivated
  // by physics analysis. Decision is put here.
  //
 
  if (isPrimary(part)) {
    return true;
  }
  if (isFromPrimaryDecayChain(part)) {
    return true;
  }
  if (isFromPrimaryPairProduction(part)) {
    return true;
  }
 
  return false;
}
 
TClonesArray* Stack::GetListOfParticles()
{
  LOG(fatal) << "Stack::GetListOfParticles interface not implemented\n";
  return nullptr;
}
 
bool Stack::isTrackDaughterOf(int trackid, int parentid) const
{
  // a daughter trackid should be larger than parentid
  if (trackid < parentid) {
    return false;
  }
 
  if (trackid == parentid) {
    return true;
  }
 
  auto mother = getMotherTrackId(trackid);
  while (mother != -1) {
    if (mother == parentid) {
      return true;
    }
    mother = getMotherTrackId(mother);
  }
  return false;
}
 
void Stack::fillParentIDs(std::vector<int>& parentids) const
{
  parentids.clear();
  int mother = mIndexOfCurrentTrack;
  do {
    if (mother != -1) {
      parentids.push_back(mother);
    }
    mother = getMotherTrackId(mother);
  } while (mother != -1);
}
 
void Stack::ReorderKine(std::vector<MCTrack>& particles, std::vector<int>& reOrderedIndices)
{
  //
  // Particles are ordered in a way that descendants of a particle appear next to each other.
  // This has the advantage that their position in the stack can be identified by two number,
  // for example the index of the first and last descentant.
  // The result of the ordering is returned via the look-up table reOrderedIndices
  //
 
  Int_t ntr = (int)(particles.size());
  std::vector<bool> done(ntr, false);
 
  int indexoffset = mTracks->size();
  Int_t index = 0;
  Int_t imoOld = 0;
  for (Int_t i = 0; i < ntr; i++) {
    reOrderedIndices[i] = i;
  }
 
  for (Int_t i = -1; i < ntr; i++) {
    if (i != -1) {
      // secondaries
      if (!done[i]) {
        reOrderedIndices[index] = i;
        index++;
        done[i] = true;
      }
      imoOld = i;
    } else {
      // current primary
      imoOld = mIndexOfCurrentPrimary - indexoffset;
    }
    for (Int_t j = i + 1; j < ntr; j++) {
      if (!done[j]) {
        if ((particles[j]).getMotherTrackId() == imoOld) {
          reOrderedIndices[index] = j;
          index++;
          done[j] = true;
        } // child found
      }   // done
    }     // j
  }       // i
}
 
FairGenericStack* Stack::CloneStack() const { return new o2::data::Stack(*this); }
 
ClassImp(o2::data::Stack);