{

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;

: 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();

LOG(info) << param;

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 (!mTransportPrimary) {

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; };

} else {

mTransportPrimary = transportPrimary;

}

: FairGenericStack(rhs),

mStack(),

mParticles(),

mTracks(nullptr),

mIndexMap(),

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),

mTrackRefs(new std::vector<o2::TrackReference>),

mIsG4Like(rhs.mIsG4Like)

{

LOG(debug) << "copy constructor called";

mTracks = new std::vector<MCTrack>();

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);

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);

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(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);

{

// 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);

}

{

// 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);

}

{

mIndexOfCurrentTrack = iTrack;

if (iTrack < mPrimaryParticles.size()) {

auto& p = mPrimaryParticles[iTrack];

mCurrentParticle = p;

mIndexOfCurrentPrimary = iTrack;

} else {

mCurrentParticle = mCurrentParticle0;

}

{

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;

{

if (mMCEventStats) {

mMCEventStats->setNHits(mHitCounter);

mMCEventStats->setNTransportedTracks(mNumberOfEntriesInParticles);

mMCEventStats->setNKeptTracks(mTracks->size());

}

{

/// This interface is not implemented since we are filtering/filling the output array

/// after each primary ... just give a summary message

LOG(info) << "Stack: " << mTracks->size() << " out of " << mNumberOfEntriesInParticles << " stored \n";

{

// 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();

{

// 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.";

{

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();

mHitCounter = 0;

{

FairRootManager::Instance()->RegisterAny("MCTrack", mTracks, kTRUE);

FairRootManager::Instance()->RegisterAny("TrackRefs", mTrackRefs, kTRUE);

{

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();

}

}

{

Int_t verbose = 0;

if (option) {

verbose = 1;

}

Print(verbose);

{

if (mIndexOfCurrentTrack < mNumberOfPrimaryParticles) {

auto& part = mTracks->at(mIndexOfCurrentTrack);

part.setHit(iDet);

} else {

Int_t iTrack = mTrackIDtoParticlesEntry[mIndexOfCurrentTrack];

auto& part = mParticles[iTrack];

part.setHit(iDet);

}

mHitCounter++;

{

mHitCounter++;

auto& part = mParticles[iTrack];

part.setHit(iDet);

{

TParticle* currentPart = GetCurrentTrack();

if (currentPart) {

return currentPart->GetFirstMother();

} else {

return -1;

}

{

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;

{

/** check if primary **/

if (part.getProcess() == kPPrimary || part.getMotherTrackId() < 0) {

return true;

}

/** not primary **/

return false;

{

/** check if the particle is from the

/** check if from decay **/

if (part.getProcess() != kPDecay) {

return false;

}

/** check if mother is primary **/

auto imother = part.getMotherTrackId();

auto mother = mParticles[imother];

if (isPrimary(mother)) {

return true;

}

/** else check if mother is from primary decay **/

return isFromPrimaryDecayChain(mother);

{

/** check if the particle is from

/** check if from pair production **/

if (part.getProcess() != kPPair) {

return false;

}

/** check if mother is primary **/

auto imother = part.getMotherTrackId();

auto mother = mParticles[imother];

if (isPrimary(mother)) {

return true;

}

/** else check if mother is from primary decay **/

return isFromPrimaryDecayChain(mother);

{

//

// 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;

{

// 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;

{

parentids.clear();

int mother = mIndexOfCurrentTrack;

do {

if (mother != -1) {

parentids.push_back(mother);

}

mother = getMotherTrackId(mother);

} while (mother != -1);

{

//

// 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