Refactor vertex selection in StPicoDstMaker

StPicoDstMaker/StPicoDstMaker.cxx

@@ -204,13 +204,10 @@ Int_t StPicoDstMaker::Init()
   {
     case PicoIoMode::IoWrite:
 
-      if (mVtxMode == PicoVtxMode::NotSet)
+      if (mVtxMode == PicoVtxMode::NotSet && setVtxModeAttr() != kStOK)
       {
-        if (setVtxModeAttr() != kStOK)
-        {
-          LOG_ERROR << "Pico Vertex Mode is not set ... " << endm;
-          return kStErr;
-        }
+        LOG_ERROR << "Pico Vertex Mode is not set ... " << endm;
+        return kStErr;
       }
 
       if (mInputFileName.Length() == 0) {
@@ -677,7 +674,7 @@ void StPicoDstMaker::fillTracks()
     }
 
     int counter = mPicoArrays[StPicoArrays::Track]->GetEntries();
-    new((*(mPicoArrays[StPicoArrays::Track]))[counter]) StPicoTrack(gTrk, pTrk, mBField, mMuDst->primaryVertex()->position(), *dcaG);
+    new((*(mPicoArrays[StPicoArrays::Track]))[counter]) StPicoTrack(gTrk, pTrk, mBField, mMuDst->primaryVertex(), *dcaG);
 
     StPicoTrack* picoTrk = (StPicoTrack*)mPicoArrays[StPicoArrays::Track]->At(counter);
 
@@ -1145,44 +1142,54 @@ bool StPicoDstMaker::selectVertex()
 {
   StMuPrimaryVertex* selectedVertex = nullptr;
 
-  if (mVtxMode == PicoVtxMode::Default)
+  switch (mVtxMode)
   {
+  case PicoVtxMode::Default:
     // choose the default vertex, i.e. the first vertex
-    mMuDst->setVertexIndex(0);
-    selectedVertex = mMuDst->primaryVertex();
-  }
-  else if (mVtxMode == PicoVtxMode::Vpd || mVtxMode == PicoVtxMode::VpdOrDefault)
-  {
-    if(mVtxMode == PicoVtxMode::VpdOrDefault)
-    {
-      mMuDst->setVertexIndex(0);
-    }
+    selectedVertex = StMuDst::setCurrentVertex(0);
+    break;
 
-    StBTofHeader const* mBTofHeader = mMuDst->btofHeader();
+  case PicoVtxMode::Vpd:
+    selectedVertex = findVpdVertex(*mMuDst);
+    StMuDst::setCurrentVertex(selectedVertex);
+    break;
 
-    if (mBTofHeader && fabs(mBTofHeader->vpdVz()) < 200)
-    {
-      float vzVpd = mBTofHeader->vpdVz();
+  case PicoVtxMode::VpdOrDefault:
+    selectedVertex = findVpdVertex(*mMuDst);
 
-      for (unsigned int iVtx = 0; iVtx < mMuDst->numberOfPrimaryVertices(); ++iVtx)
-      {
-        StMuPrimaryVertex* vtx = mMuDst->primaryVertex(iVtx);
-        if (!vtx) continue;
+    if (selectedVertex)
+      StMuDst::setCurrentVertex(selectedVertex);
+    else
+      selectedVertex = StMuDst::setCurrentVertex(0);
 
-        if (fabs(vzVpd - vtx->position().z()) < 3.)
-        {
-          mMuDst->setVertexIndex(iVtx);
-          selectedVertex = mMuDst->primaryVertex();
-          break;
-        }
-      }
-    }
-  }
-  else // default case
-  {
+    break;
+
+  default:
     LOG_ERROR << "Pico Vtx Mode not set!" << endm;
   }
 
   // Retrun false if selected vertex is not valid
   return selectedVertex ? true : false;
 }
+
+
+StMuPrimaryVertex* StPicoDstMaker::findVpdVertex(const StMuDst& muDst) const
+{
+  StBTofHeader const* mBTofHeader = muDst.btofHeader();
+
+  if (!mBTofHeader || fabs(mBTofHeader->vpdVz()) >= 200)
+    return nullptr;
+
+  float vzVpd = mBTofHeader->vpdVz();
+
+  for (unsigned int iVtx = 0; iVtx < muDst.numberOfPrimaryVertices(); ++iVtx)
+  {
+    StMuPrimaryVertex* vtx = muDst.primaryVertex(iVtx);
+
+    if (!vtx || fabs(vzVpd - vtx->position().z()) >= 3.) continue;
+
+    return vtx;
+  }
+
+  return nullptr;
+}

StPicoDstMaker/StPicoDstMaker.h

@@ -15,6 +15,7 @@ class TFile;
 class TTree;
 class StMuDst;
 class StMuTrack;
+class StMuPrimaryVertex;
 class StEmcCollection;
 class StEmcPosition;
 class StEmcGeom;
@@ -28,7 +29,21 @@ class StPicoDstMaker : public StMaker
 {
 public:
   enum PicoIoMode {IoWrite=1, IoRead=2};
-  enum PicoVtxMode {NotSet=0, Default=1, Vpd=2, VpdOrDefault=3};
+
+  /// Various selection criteria for picking a vertex from the list of available
+  /// vertices in muDst. If requested vertex not found the event is skipped
+  enum PicoVtxMode
+  {
+    /// Default "unspecified" mode to force user pick one of the modes below
+    NotSet = 0,
+    /// Select the first vertex in the muDst vertex collection
+    Default = 1,
+    /// Select first available vertex within a window around the VPD one
+    Vpd = 2,
+    /// If a vertex matching the VPD one is not found pick the first one from
+    /// the muDst collection
+    VpdOrDefault = 3
+  };
 
   StPicoDstMaker(char const* name = "PicoDst");
   StPicoDstMaker(PicoIoMode ioMode, char const* fileName = "", char const* name = "PicoDst");
@@ -128,6 +143,10 @@ class StPicoDstMaker : public StMaker
   /// A pointer to the main input/outpur picoDst structure containing all `TObjArray`s
   StPicoDst*  mPicoDst;
 
+  /// Finds a primary vertex in the list of muDst vertices that matches the VPD
+  /// vertex
+  StMuPrimaryVertex* findVpdVertex(const StMuDst& muDst) const;
+
   StEmcCollection* mEmcCollection;
   StEmcPosition*   mEmcPosition;
   StEmcGeom*       mEmcGeom[4];

StPicoEvent/StPicoEvent.cxx

@@ -42,11 +42,11 @@ StPicoEvent::StPicoEvent(StMuDst const& muDst) : StPicoEvent()
   mFillId = ev->runInfo().beamFillNumber(blue);
   mBField = ev->magneticField();
 
-  mPrimaryVertex = ev->primaryVertexPosition();
-  mPrimaryVertexError = ev->primaryVertexErrors();
-
   if (StMuPrimaryVertex* pv = muDst.primaryVertex())
   {
+    mPrimaryVertex = ev->primaryVertexPosition();
+    mPrimaryVertexError = ev->primaryVertexErrors();
+
     mRanking = pv->ranking();
     mNBEMCMatch = pv->nBEMCMatch();
     mNBTOFMatch = pv->nBTOFMatch();

StPicoEvent/StPicoTrack.cxx

@@ -6,6 +6,7 @@
 #include "StPicoEvent/StPicoTrack.h"
 #include "St_base/StMessMgr.h"
 #include "StMuDSTMaker/COMMON/StMuTrack.h"
+#include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
 
 
 //----------------------------------------------------------------------------------
@@ -23,7 +24,7 @@ StPicoTrack::StPicoTrack() : TObject(),
 }
 
 //----------------------------------------------------------------------------------
-StPicoTrack::StPicoTrack(StMuTrack const* const gTrk, StMuTrack const* const pTrk, double const B, StThreeVectorD const& pVtx, StDcaGeometry const& dcaG)
+StPicoTrack::StPicoTrack(StMuTrack const* const gTrk, StMuTrack const* const pTrk, double const B, const StMuPrimaryVertex* pVtx, StDcaGeometry const& dcaG)
   : StPicoTrack()
 {
   if (!gTrk || gTrk->type() != global || (pTrk && (pTrk->type() != primary || pTrk->id() != gTrk->id())))
@@ -41,11 +42,15 @@ StPicoTrack::StPicoTrack(StMuTrack const* const gTrk, StMuTrack const* const pTr
     mPMomentum = pTrk->p();
   }
 
-  // Calculate global momentum and position at point of DCA to the pVtx
-  StPhysicalHelixD gHelix = dcaG.helix();
-  gHelix.moveOrigin(gHelix.pathLength(pVtx));
-  mGMomentum = gHelix.momentum(B * kilogauss);
-  mOrigin = gHelix.origin();
+  // Calculate global momentum and position at point of DCA to the pVtx if the
+  // vertex is valid
+  if (pVtx)
+  {
+    StPhysicalHelixD gHelix = dcaG.helix();
+    gHelix.moveOrigin(gHelix.pathLength( pVtx->position() ));
+    mGMomentum = gHelix.momentum(B * kilogauss);
+    mOrigin = gHelix.origin();
+  }
 
   mDedx      = gTrk->dEdx() * 1.e6;
   mDnDx      = gTrk->probPidTraits().dNdxFit();

StPicoEvent/StPicoTrack.h

@@ -9,6 +9,8 @@
 #include "StarClassLibrary/SystemOfUnits.h"
 
 
+class StDcaGeometry;
+class StMuPrimaryVertex;
 class StMuTrack;
 class StDcaGeometry;
 
@@ -21,7 +23,7 @@ class StPicoTrack : public TObject
 
   /// ctor. Note: primary track should be associated with the StPicoEvent::mPrimaryVertex
   StPicoTrack(StMuTrack const* globalTrack, StMuTrack const* primaryTrack,
-     double magField, StThreeVectorD const& pVtx, StDcaGeometry const& dcaG);
+     double magField, const StMuPrimaryVertex* pVtx, StDcaGeometry const& dcaG);
 
   virtual ~StPicoTrack() {}