Revert "Revert "Implementation of simulation inhibit-pause-trigger-continue features (#5305)""

This reverts commit 27ac59a.
It re-enables the inhibit-pause-trigger-continue feature
and fixes a problem in initialization order of stack and parameters.

Author: sawenzel

DataFormats/simulation/include/SimulationDataFormat/MCTrack.h

@@ -329,6 +329,8 @@ inline MCTrackT<T>::MCTrackT(const TParticle& part)
 {
   // our convention is to communicate the process as (part) of the unique ID
   setProcess(part.GetUniqueID());
+  // extract storage flag
+  setStore(part.TestBit(ParticleStatus::kKeep));
   // extract toBeDone flag
   setToBeDone(part.TestBit(ParticleStatus::kToBeDone));
   // extract inhibited flag

DataFormats/simulation/include/SimulationDataFormat/Stack.h

@@ -90,7 +90,8 @@ class Stack : public FairGenericStack
 
   void 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 proc, Int_t& ntr, Double_t weight, Int_t is, Int_t secondParentId, Int_t daughter1Id, Int_t daughter2Id,
+                 TMCProcess proc2);
 
   // similar function taking a particle
   void PushTrack(Int_t toBeDone, TParticle&);

DataFormats/simulation/src/Stack.cxx

@@ -103,6 +103,7 @@ Stack::Stack(Int_t size)
     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;
   }
@@ -183,12 +184,13 @@ void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px
                       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);
+  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 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: ",
@@ -214,36 +216,38 @@ void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px
   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 (proc == kPPrimary) {
+  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);
-    p.SetBit(ParticleStatus::kPrimary);
-    if (toBeDone == 1) {
-      handleTransportPrimary(p);
-    } else {
-      p.SetBit(ParticleStatus::kToBeDone, 0);
+    p.SetBit(ParticleStatus::kKeep, 1);
+    if (p.TestBit(ParticleStatus::kToBeDone)) {
+      mNumberOfPrimariesforTracking++;
     }
     mNumberOfPrimaryParticles++;
     mPrimaryParticles.push_back(p);
     mTracks->emplace_back(p);
   } else {
-    p.SetBit(ParticleStatus::kPrimary, 0);
-    if (toBeDone == 1) {
-      p.SetBit(ParticleStatus::kToBeDone, 1);
-    } else {
-      p.SetBit(ParticleStatus::kToBeDone, 0);
-    }
     mParticles.emplace_back(p);
     mCurrentParticle0 = p;
   }
@@ -255,10 +259,11 @@ void Stack::handleTransportPrimary(TParticle& p)
   // this function tests whether we really want to transport
   // this particle and sets the relevant bits accordingly
 
-  if (mTransportPrimary(p, mPrimaryParticles)) {
-    p.SetBit(ParticleStatus::kToBeDone, 1);
-    mNumberOfPrimariesforTracking++;
-  } else {
+  if (!p.TestBit(ParticleStatus::kToBeDone) || !p.TestBit(ParticleStatus::kPrimary)) {
+    return;
+  }
+
+  if (!mTransportPrimary(p, mPrimaryParticles)) {
     p.SetBit(ParticleStatus::kToBeDone, 0);
     p.SetBit(ParticleStatus::kInhibited, 1);
   }
@@ -271,15 +276,15 @@ void Stack::PushTrack(int toBeDone, TParticle& p)
   // This method is called
   //
   // - during parallel simulation to push primary particles (called by the stack itself)
-  if (p.GetUniqueID() == 0) {
+  if (p.TestBit(ParticleStatus::kPrimary)) {
     // one to one mapping for primaries
     mIndexMap[mNumberOfPrimaryParticles] = mNumberOfPrimaryParticles;
-    // Push particle on the stack
-    if (p.TestBit(ParticleStatus::kPrimary) && p.TestBit(ParticleStatus::kToBeDone)) {
-      handleTransportPrimary(p);
-    }
     mNumberOfPrimaryParticles++;
     mPrimaryParticles.push_back(p);
+    // Push particle on the stack
+    if (p.TestBit(ParticleStatus::kToBeDone)) {
+      mNumberOfPrimariesforTracking++;
+    }
     mStack.push(p);
     mTracks->emplace_back(p);
   }

Generators/CMakeLists.txt

@@ -27,6 +27,7 @@ o2_add_library(Generators
                        src/GeneratorTGenerator.cxx
                        src/GeneratorExternalParam.cxx
                        src/GeneratorFromFile.cxx
+                       src/GeneratorFromO2KineParam.cxx
                        src/PDG.cxx
                        src/PrimaryGenerator.cxx
                        src/InteractionDiamondParam.cxx
@@ -68,6 +69,7 @@ set(headers
     include/Generators/GeneratorTGenerator.h
     include/Generators/GeneratorExternalParam.h
     include/Generators/GeneratorFromFile.h
+    include/Generators/GeneratorFromO2KineParam.h
     include/Generators/PDG.h
     include/Generators/PrimaryGenerator.h
     include/Generators/InteractionDiamondParam.h

Generators/include/Generators/GeneratorFromFile.h

@@ -78,11 +78,15 @@ class GeneratorFromO2Kine : public o2::eventgen::Generator
   void SetStartEvent(int start);
 
  private:
+  /** methods that can be overridden **/
+  void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override;
+
   TFile* mEventFile = nullptr;     //! the file containing the persistent events
   TBranch* mEventBranch = nullptr; //! the branch containing the persistent events
   int mEventCounter = 0;
   int mEventsAvailable = 0;
   bool mSkipNonTrackable = true; //! whether to pass non-trackable (decayed particles) to the MC stack
+  bool mContinueMode = false;    //! whether we want to continue simulation of previously inhibited tracks
   ClassDefOverride(GeneratorFromO2Kine, 1);
 };
 

Generators/include/Generators/GeneratorFromO2KineParam.h

@@ -0,0 +1,39 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// 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.
+
+/// \author R+Preghenella - January 2021
+
+#ifndef ALICEO2_EVENTGEN_GENERATORFROMO2KINEPARAM_H_
+#define ALICEO2_EVENTGEN_GENERATORFROMO2KINEPARAM_H_
+
+#include "CommonUtils/ConfigurableParam.h"
+#include "CommonUtils/ConfigurableParamHelper.h"
+
+namespace o2
+{
+namespace eventgen
+{
+
+/**
+ ** a parameter class/struct to keep the settings of
+ ** the FromO2Kine event generator and
+ ** allow the user to modify them 
+ **/
+
+struct GeneratorFromO2KineParam : public o2::conf::ConfigurableParamHelper<GeneratorFromO2KineParam> {
+  bool skipNonTrackable = true;
+  bool continueMode = false;
+  O2ParamDef(GeneratorFromO2KineParam, "GeneratorFromO2Kine");
+};
+
+} // end namespace eventgen
+} // end namespace o2
+
+#endif // ALICEO2_EVENTGEN_GENERATORFROMO2KINEPARAM_H_

Generators/src/Generator.cxx

@@ -131,7 +131,8 @@ Bool_t
                         particle.GetStatusCode() == 1,
                         particle.Energy() * mEnergyUnit,
                         particle.T() * mTimeUnit,
-                        particle.GetWeight());
+                        particle.GetWeight(),
+                        (TMCProcess)particle.GetUniqueID());
   }
 
   /** success **/

Generators/src/GeneratorFromFile.cxx

@@ -9,7 +9,9 @@
 // or submit itself to any jurisdiction.
 
 #include "Generators/GeneratorFromFile.h"
+#include "Generators/GeneratorFromO2KineParam.h"
 #include "SimulationDataFormat/MCTrack.h"
+#include "SimulationDataFormat/MCEventHeader.h"
 #include <FairLogger.h>
 #include <FairPrimaryGenerator.h>
 #include <TBranch.h>
@@ -176,6 +178,13 @@ GeneratorFromO2Kine::GeneratorFromO2Kine(const char* name)
     }
   }
   LOG(ERROR) << "Problem reading events from file " << name;
+
+  // read and set params
+  auto& param = GeneratorFromO2KineParam::Instance();
+  LOG(INFO) << "Instance \'FromO2Kine\' generator with following parameters";
+  LOG(INFO) << param;
+  mSkipNonTrackable = param.skipNonTrackable;
+  mContinueMode = param.continueMode;
 }
 
 void GeneratorFromO2Kine::SetStartEvent(int start)
@@ -201,8 +210,9 @@ bool GeneratorFromO2Kine::importParticles()
     mEventBranch->GetEntry(mEventCounter);
 
     for (auto& t : *tracks) {
-      // I guess we only want primaries (unless later on we continue a simulation)
-      if (!t.isPrimary()) {
+
+      // in case we do not want to continue, take only primaries
+      if (!mContinueMode && !t.isPrimary()) {
         continue;
       }
 
@@ -221,9 +231,16 @@ bool GeneratorFromO2Kine::importParticles()
       auto vt = t.T();
       auto weight = 1.; // p.GetWeight() ??
       auto wanttracking = t.getToBeDone();
+
+      if (mContinueMode) { // in case we want to continue, do only inhibited tracks
+        wanttracking &= t.getInhibited();
+      }
+
       LOG(DEBUG) << "Putting primary " << pdg;
 
       mParticles.push_back(TParticle(pdg, wanttracking, m1, m2, d1, d2, px, py, pz, e, vx, vy, vz, vt));
+      mParticles.back().SetUniqueID((unsigned int)t.getProcess()); // we should propagate the process ID
+
       particlecounter++;
     }
     mEventCounter++;
@@ -240,6 +257,17 @@ bool GeneratorFromO2Kine::importParticles()
   return false;
 }
 
+void GeneratorFromO2Kine::updateHeader(o2::dataformats::MCEventHeader* eventHeader)
+{
+  /** update header **/
+
+  // put information about input file and event number of the current event
+
+  eventHeader->putInfo<std::string>("generator", "generatorFromO2Kine");
+  eventHeader->putInfo<std::string>("inputFile", mEventFile->GetName());
+  eventHeader->putInfo<int>("inputEventNumber", mEventCounter - 1);
+}
+
 } // namespace eventgen
 } // end namespace o2
 

Generators/src/GeneratorFromO2KineParam.cxx

@@ -0,0 +1,14 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// 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.
+
+/// \author R+Preghenella - January 2021
+
+#include "Generators/GeneratorFromO2KineParam.h"
+O2ParamImpl(o2::eventgen::GeneratorFromO2KineParam);

Generators/src/GeneratorsLinkDef.h

@@ -48,6 +48,8 @@
 #endif
 #pragma link C++ class o2::eventgen::GeneratorFromFile + ;
 #pragma link C++ class o2::eventgen::GeneratorFromO2Kine + ;
+#pragma link C++ class o2::eventgen::GeneratorFromO2KineParam + ;
+#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GeneratorFromO2KineParam> + ;
 #pragma link C++ class o2::PDG + ;
 #pragma link C++ class o2::eventgen::PrimaryGenerator + ;