Merge branch 'dev' of https://github.com/ExpertRootGroup/er into dev

Gadast/ERGadast.cxx

@@ -109,7 +109,7 @@ int ERGadast::ParentGammaTrackId(int track_id) const {
 
   auto* track = track_it->second;
   auto* parent_track = track;
-  while (parent_track->GetPdgCode() != 22 || parent_track->GetMotherId() != -1) {
+  while (parent_track->GetPdgCode() != 22 && parent_track->GetMotherId() != -1) {
     parent_track = id_to_track[parent_track->GetMotherId()];
   }
 

Gadast/ERGadast.h

@@ -147,7 +147,7 @@ class ERGadast : public ERDetector
   TClonesArray*  fCsIPoints = nullptr;                  ///< The CsI point collection
   TClonesArray*  fLaBrPoints = nullptr;                 ///< The LaBr point collection
   TClonesArray*  fGadastSteps = nullptr;                ///< The all steps collection
-  TClonesArray*  fMCTracks = nullptr;                   ///< The all steps collection
+  TClonesArray*  fMCTracks = nullptr;                   ///< The mc tracks collection
 
   Bool_t fStoreSteps;                         ///< Flag to storing all steps in sensetive volume
 

ND/CMakeLists.txt

@@ -12,6 +12,7 @@ ${CMAKE_SOURCE_DIR}/BeamDet/data
 ${CMAKE_SOURCE_DIR}/data/
 ${CMAKE_SOURCE_DIR}/ND/
 ${CMAKE_SOURCE_DIR}/ND/data
+${CMAKE_SOURCE_DIR}/decayers
 )
 
 include_directories( ${INCLUDE_DIRECTORIES})

ND/ERND.cxx

@@ -8,13 +8,21 @@
 
 #include "ERND.h"
 
+#include <unordered_map>
+
 #include "TClonesArray.h"
 #include "TParticle.h"
 #include "TVirtualMC.h"
 #include "TString.h"
 
 #include "FairRootManager.h"
 #include "FairLogger.h"
+#include "FairRunSim.h"
+
+#include "ERMCTrack.h"
+#include "ERStack.h"
+#include "ERDecayMCEventHeader.h"
+
 //--------------------------------------------------------------------------------------------------
 ERND::ERND()
   : FairDetector("ERND", kTRUE)
@@ -76,6 +84,9 @@ Bool_t ERND::ProcessHits(FairVolume* vol) {
   static Int_t          stilbenNr;
   static Double_t       lightYield;
 
+  static Int_t          parentTrackId = -1;
+  static Int_t          parentPdg = -1;
+
   gMC->SetMaxStep(fStep);
 
   if ( gMC->IsTrackEntering() ) { // Return true if this is the first step of the track in the current volume
@@ -118,14 +129,17 @@ Bool_t ERND::ProcessHits(FairVolume* vol) {
 	{ 
     gMC->TrackPosition(posOut);
     gMC->TrackMomentum(momOut);
-
-	  if (eLoss > 0. && gMC->TrackCharge()!=0){
+
+	  if (eLoss > 0. && gMC->TrackCharge()!=0)
+    {
+      FindParentParticle(trackID, parentTrackId, parentPdg);
       AddPoint( eventID, trackID, mot0TrackID, pdg,
                 TVector3(posIn.X(),   posIn.Y(),   posIn.Z()),
                 TVector3(posOut.X(),  posOut.Y(),  posOut.Z()),
                 TVector3(momIn.Px(),  momIn.Py(),  momIn.Pz()),
                 TVector3(momOut.Px(), momOut.Py(), momOut.Pz()),
-                time, length, eLoss, stilbenNr, lightYield);
+                time, length, eLoss, stilbenNr, lightYield,
+                parentTrackId, parentPdg);
     }
   }
   return kTRUE;
@@ -140,7 +154,8 @@ void ERND::EndOfEvent() {
 void ERND::Register() {
   FairRootManager* ioman = FairRootManager::Instance();
   if (!ioman)
-	Fatal("Init", "IO manager is not set");	
+	Fatal("Init", "IO manager is not set");
+  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
   ioman->Register("NDPoint","ND", fNDPoints, kTRUE);
 }
 //--------------------------------------------------------------------------------------------------
@@ -178,18 +193,69 @@ void ERND::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) {
   LOG(DEBUG) << "decector: " << cl2->GetEntriesFast() << " merged entries" << FairLogger::endl;
 }
 //--------------------------------------------------------------------------------------------------
-ERNDPoint* ERND::AddPoint(Int_t eventID, Int_t trackID,
-				    Int_t mot0trackID,
-				    Int_t pdg,
-				    TVector3 posIn,
-				    TVector3 posOut, TVector3 momIn,
-				    TVector3 momOut, Double_t time,
-				    Double_t length, Double_t eLoss, Int_t stilbenNr, Float_t lightYield) {
+ERNDPoint* ERND::AddPoint(
+    Int_t eventID, Int_t trackID, Int_t mot0trackID, Int_t pdg,
+    TVector3 posIn, TVector3 posOut, TVector3 momIn,
+    TVector3 momOut, Double_t time, Double_t length,
+    Double_t eLoss, Int_t stilbenNr, Float_t lightYield,
+    Int_t parentTrackId, Int_t parentPdg) 
+{
   TClonesArray& clref = *fNDPoints;
   Int_t size = clref.GetEntriesFast();
   return new(clref[size]) ERNDPoint(eventID, trackID, mot0trackID, pdg,
-					  posIn, posOut, momIn, momOut, time, length, eLoss,stilbenNr,lightYield);
-
+					  posIn, posOut, momIn, momOut, time, length, eLoss,stilbenNr, lightYield,
+            parentTrackId, parentPdg);
+}
+//--------------------------------------------------------------------------------------------------
+void ERND::FindParentParticle(const int track_id, int& parentTrackId, int& parentPdg)
+{
+  dynamic_cast<ERStack*>(gMC->GetStack())->FillTrackArray();
+  std::unordered_map<int, ERMCTrack*> id_to_track;
+  for (int i(0); i < fMCTracks->GetEntriesFast(); i++) {
+    auto* track = dynamic_cast<ERMCTrack*>(fMCTracks->At(i));
+    id_to_track[track->Id()] = track;
+  }
+
+  auto track_it = id_to_track.find(track_id);
+  if (track_it == id_to_track.end())
+    LOG(FATAL) << "ND: Track not found in stack!" << FairLogger::endl;
+
+  auto* track = track_it->second;
+  auto* parent_track = track;
+
+  FairRunSim* run = FairRunSim::Instance();
+  if (auto* header = dynamic_cast<ERDecayMCEventHeader*>(run->GetMCEventHeader()))
+  {
+    TArrayI reaction_track_ids = header->GetOutputTrackID();
+    const auto found_in_reaction = [&reaction_track_ids](const int id)
+    {
+      for (int i = 0; i < reaction_track_ids.GetSize(); ++i)
+      {
+        if (reaction_track_ids[i] == id)
+          return true;
+      }
+      return false;
+    };
+
+    while (!found_in_reaction(parent_track->Id()) && parent_track->GetMotherId() != -1)
+    {
+      parent_track = id_to_track[parent_track->GetMotherId()];
+    } 
+  }
+  else
+  {
+    while ((fCandidatesForParentPdgs.find(parent_track->GetPdgCode()) != fCandidatesForParentPdgs.end())
+        && parent_track->GetMotherId() != -1)
+    {
+      parent_track = id_to_track[parent_track->GetMotherId()];
+    }
+  }
+
+  parentTrackId = parent_track->Id();
+  if (parent_track->GetMotherId() != -1)
+  {
+    parentPdg = parent_track->GetPdgCode();
+  }
 }
 //--------------------------------------------------------------------------------------------------
 void ERND::ConstructGeometry() {

ND/ERND.h

@@ -9,6 +9,8 @@
 #ifndef ERND_H
 #define ERND_H
 
+#include <set>
+
 #include "ERDetector.h"
 #include "ERNDPoint.h"
 
@@ -33,6 +35,9 @@ class ERND : public FairDetector {
   /** @brief Set propagation step in sensetive volumes.
    ** Special process stepLimiter must be in TG4RunConfiguration **/
   void SetMaxStep(Double_t step) {fStep = step;}
+
+  void SetCandidatesForParentPDG(const std::set<Int_t>& pdgs) {fCandidatesForParentPdgs = pdgs;}
+
   /** Virtual method ProcessHits
    ** Defines the action to be taken when a step is inside the
    ** active volume. Creates a ERNDPoint and adds it to the
@@ -68,16 +73,23 @@ class ERND : public FairDetector {
   virtual Bool_t CheckIfSensitive(std::string name);
 private:
   TClonesArray* fNDPoints = nullptr; //!  The point collection
+  TClonesArray*  fMCTracks = nullptr;
   Double_t fStep = 1.; //! Max lengt of step of track propagetion in sensetive volume
+  std::set<Int_t> fCandidatesForParentPdgs = {22, 2112};
 private:
+  void FindParentParticle(int track_id, int& parentTrackId, int& parentPdg);
   /// Adds a NeuRadPoint to the Point Collection
-  ERNDPoint* AddPoint(Int_t eventID, Int_t trackID,
-			  Int_t mot0trackID,
-			  Int_t pdg,
-			  TVector3 posIn,
-			  TVector3 pos_out, TVector3 momIn,
-			  TVector3 momOut, Double_t time,
-			  Double_t length, Double_t eLoss, Int_t stilbenNr,Float_t lightYield);
+  ERNDPoint* AddPoint(
+    Int_t eventID, Int_t trackID,
+    Int_t mot0trackID,
+    Int_t pdg,
+    TVector3 posIn,
+    TVector3 pos_out, TVector3 momIn,
+    TVector3 momOut, Double_t time,
+    Double_t length, Double_t eLoss, 
+    Int_t stilbenNr,Float_t lightYield,
+    Int_t parentTrackId, Int_t parentPdg
+  );
   ClassDef(ERND,1);
 };
 

ND/data/ERNDPoint.cxx

@@ -24,13 +24,16 @@ ERNDPoint::ERNDPoint(Int_t eventID, Int_t trackID,
 		  Int_t pdg,
 		  TVector3 posIn,
 		  TVector3 posOut, TVector3 momIn, TVector3 momOut,
-		  Double_t tof, Double_t length, Double_t eLoss, Int_t stilbenNr,Float_t lightYield)
+		  Double_t tof, Double_t length, Double_t eLoss,
+      Int_t stilbenNr,Float_t lightYield,
+      Int_t parentTrackId, Int_t parentPdgId)
   : FairMCPoint(trackID, -1., posIn, momIn, tof, length, eLoss),
     fMot0TrackID(mot0trackID),
     fEventID(eventID),
     fX_out(posOut.X()), fY_out(posOut.Y()), fZ_out(posOut.Z()),
     fPx_out(momOut.X()), fPy_out(momOut.Y()), fPz_out(momOut.Z()),
-    fStilbenNr(stilbenNr),fLightYield(lightYield), fPdg(pdg)
+    fStilbenNr(stilbenNr),fLightYield(lightYield), fPdg(pdg),
+    fParentTrackId(parentTrackId), fParentPdgId(parentPdgId)
 {
 }
 // -------------------------------------------------------------------------

ND/data/ERNDPoint.h

@@ -39,12 +39,23 @@ class ERNDPoint : public FairMCPoint
    *@param length   Track length since creation [cm]
    *@param eLoss    Energy deposit [MeV]
    **/
-  ERNDPoint(Int_t eventID, Int_t trackID,
-		  Int_t mot0trackID,
-		  Int_t pdg,
-		  TVector3 posIn,
-		  TVector3 posOut, TVector3 momIn, TVector3 momOut,
-		  Double_t tof, Double_t length, Double_t eLoss, Int_t stilbenNr, Float_t lightYield);
+  ERNDPoint(
+    Int_t eventID,
+    Int_t trackID,
+    Int_t mot0trackID,
+    Int_t pdg,
+    TVector3 posIn,
+    TVector3 posOut,
+    TVector3 momIn,
+    TVector3 momOut,
+    Double_t tof,
+    Double_t length,
+    Double_t eLoss,
+    Int_t stilbenNr,
+    Float_t lightYield,
+    Int_t parentTrackId,
+    Int_t parentPdgId
+  );
 
 
   /** Copy constructor **/
@@ -101,7 +112,9 @@ class ERNDPoint : public FairMCPoint
   Int_t fStilbenNr;
   Float_t fLightYield;
   Int_t fPdg;
-
+  Int_t fParentTrackId = -1;                      ///< id of parrent particle track.
+  Int_t fParentPdgId = -1;                        ///< pdg of parrent particle.
+
   ClassDef(ERNDPoint,1)
 };
 #endif

README.md

@@ -37,7 +37,7 @@ for different studies.
 Due to high dependencies on the external packages of various versions, we strongly recommend to use
 the docker-image [docker](https://www.docker.com) to deploy ER installation.
 
-1. Install docker desktop on your system: https://www.docker.com/get-started
+1. Install docker engine on your system: https://docs.docker.com/engine/install/ubuntu/
 2. Follow docker post-installation steps: https://docs.docker.com/engine/install/linux-postinstall/
 2. Clone ER repository:
 
@@ -52,36 +52,13 @@ git checkout dev
 docker build --build-arg ER=dev -t er .
 ```
 
-4. Enjoy  _ROOT_, _Go4_, _AccDAQ_ and _ER_:
+4. Run _er_ container, compile updated sources, run simulation:
 
-* To run simulation in interactive session:
-
-```
-docker run -v /home/vitaliy/er/macro/EXP1904_H7:/opt/run -w /opt/run -e DISPLAY=$DISPLAY -it er /bin/bash
-root -l -q create_passive_component.C
-root -l -q create_target_exp1904.C
-root -l -q sim_digi.C
-root -l sim_digi.root
-```
-where -v is used to map your host working directory '/home/vitaliy/er/macro/EXP1904_H7' and working 
-directory  '/opt/run' inside container; -w /opt/run - set working directory for interactive session;
--e DISPLAY=$DISPLAY is needed to forward gui from container to host machine; -it - to set interactive session.
-
-* To run macro in batch mode:
-
-```
-docker run -v /Users/vitaliy/er/macro/EXP1904_H7:/opt/run -w /opt/run  er root -l -b -q create_passive_component.C
 ```
-
-5. _ER_ developing:
-
-Map _ER_ source directory on host machine and container to compile updated sources:
-
-```
-docker run --entrypoint /bin/bash -v /home/vitaliy/er:/opt/er -v /home/vitaliy/er/macro/EXP1904_H7:/opt/run -w /opt/run -e DISPLAY=$DISPLAY -it er
+docker run --entrypoint /bin/bash --net=host -v /home/vitaliy/er:/opt/er -v /home/vitaliy/er/macro/EXP1904_H7:/opt/run -w /opt/run -v /tmp/.X11-unix:/tmp/.X11-unix  -v $HOME/.Xauthority:/home/jovyan/.Xauthority:rw -e DISPLAY=$DISPLAY -it er
 cd /opt/er
 mkdir build
-cd /build
+cd build
 export SIMPATH=/opt/FairSoft/
 export FAIRROOTPATH=/opt/FairRoot/
 cmake ../ -DACCULINNA_GO4=/opt/accdaq/install/
@@ -92,6 +69,10 @@ cd /opt/run/
 root -l sim_digi.C
 ```
 
+where -v is used to map your host working directory '/home/vitaliy/er/macro/EXP1904_H7' and working 
+directory  '/opt/run' inside container; -w /opt/run - set working directory for interactive session;
+-e DISPLAY=$DISPLAY is needed to forward gui from container to host machine; -it - to set interactive session.
+
 Redifinition of `entrypoint` is needed because `er/build/config.sh` may not exists in mapped volume.
 
 ## Step by Step installation

data/ERDecayMCEventHeader.h

@@ -18,8 +18,8 @@
 #include "FairMCEventHeader.h"
 
 class ERDecayMCEventHeader : public FairMCEventHeader {
-private:
-  TVector3 fDecayPos;
+protected:
+  TVector3 fReactionPos;
   Int_t fInputTrackID;
   TArrayI fOutputTrackID;
   Int_t fOutputsTracksNb = 0;
@@ -28,7 +28,7 @@ class ERDecayMCEventHeader : public FairMCEventHeader {
 public:
   ERDecayMCEventHeader();
 
-  void SetDecayPos(TVector3 pos){fDecayPos = pos;}
+  void SetReactionPos(TVector3 pos){fReactionPos = pos;}
   void SetInputIon(Int_t iID){fInputTrackID = iID;}
   void SetTrigger(Int_t trigger){fTrigger = trigger;}
   void SetTriggerPriority(Int_t priority){fTriggerPriority = priority;}
@@ -38,7 +38,7 @@ class ERDecayMCEventHeader : public FairMCEventHeader {
 
   Int_t GetTrigger() const {return fTrigger;}
   Int_t GetTriggerPriority() const {return fTriggerPriority;}
-
+  TArrayI GetOutputTrackID() const {return fOutputTrackID;}
   void Clear();
 
   ClassDef(ERDecayMCEventHeader,1)

decayers/CMakeLists.txt

@@ -44,6 +44,6 @@ CHANGE_FILE_EXTENSION(*.cxx *.h HEADERS "${SRCS}")
 
 Set(LINKDEF ERDecayersLinkDef.h)
 Set(LIBRARY_NAME ERDecayers)
-Set(DEPENDENCIES Base Core)
+Set(DEPENDENCIES Base Core ERData)
 
 GENERATE_LIBRARY()

decayers/ER10Heto8HeEventHeader.cxx

@@ -6,10 +6,11 @@
 #include "FairLogger.h"
 
 void ER10Heto8HeEventHeader::SetData(const TVector3& position, const TLorentzVector& He8b,
-                                   const TLorentzVector& H3, const TLorentzVector& H1,
-                                   const TLorentzVector& He8d,
-                                   const TLorentzVector& He10, const TLorentzVector& n1,
-                                   const TLorentzVector& n2, const float time) {
+                                    const TLorentzVector& H3, const TLorentzVector& H1,
+                                    const TLorentzVector& He8d,
+                                    const TLorentzVector& He10, const TLorentzVector& n1,
+                                    const TLorentzVector& n2, const float time) 
+{
   fReactionPos = position;
   fHe8b = He8b;
   fH3 = H3;
@@ -22,6 +23,7 @@ void ER10Heto8HeEventHeader::SetData(const TVector3& position, const TLorentzVec
 }
 // -------------------------------------------------------------------------
 void ER10Heto8HeEventHeader::Clear() {
+  ERDecayMCEventHeader::Clear();
   fHe8b.SetXYZM(0, 0, 0, 0);
   fH3.SetXYZM(0, 0, 0, 0);
   fH1.SetXYZM(0, 0, 0, 0);