Fwd picodst integration

StRoot/StBFChain/BigFullChain.h

@@ -1949,6 +1949,7 @@ Bfc_st BFC[] = { // standard chains
   {"PicoVtxVpdOrDefault","","","-PicoVtxDefault"    ,"","","pico Vtx cut on Tof and VPD or default",kFALSE},
   {"PicoVtxFXT"     ,"","","-PicoVtxDefault"    ,"" ,"","pico Vtx constraint on FXT [198,202] mode",kFALSE},
   {"PicoVtxMtd"     ,"","","-PicoVtxDefault"             ,"" ,"","pico Vtx using MTD matching mode",kFALSE},
+  {"PicoVtxless"    ,"","","-PicoVtxDefault"             ,"" ,"","pico Vtx NOT required for FWD"   ,kFALSE},
   {"PicoCovMtxSkip" ,"","",""       ,"" ,"","Do not write covariance matrices to picoDst (default)",kFALSE},
   {"PicoCovMtxWrite","","","-PicoCovMtxSkip"   ,"" ,"","Write track covariance matrices to picoDst",kFALSE},
   {"PicoBEmcSmdSkip" ,"","",""                     ,"" ,"","Do not write BSMD to picoDst (default)",kFALSE},

StRoot/StFwdTrackMaker/macro/viz2.C

@@ -1,4 +1,4 @@
-//usr/bin/env root -l -b -q  $0'('$1')'; exit $?
+//usr/bin/env root -l -b -q  $0'("'$1'")'; exit $?
 #include "TFile.h"
 #include "TH1F.h"
 #include "TH2F.h"

StRoot/StPicoDstMaker/StPicoDstMaker.cxx

@@ -40,6 +40,10 @@
 #include "StEvent/StTriggerData.h"
 #include "StEvent/StEnumerations.h"
 #include "StEvent/StL0Trigger.h"
+#include "StEvent/StFwdTrackCollection.h"
+#include "StEvent/StFwdTrack.h"
+#include "StEvent/StFcsCluster.h"
+#include "StEvent/StPrimaryVertex.h"
 
 #include "StMuDSTMaker/COMMON/StMuDst.h"
 #include "StMuDSTMaker/COMMON/StMuEvent.h"
@@ -59,6 +63,9 @@
 #include "StMuDSTMaker/COMMON/StMuETofHeader.h"
 #include "StMuDSTMaker/COMMON/StMuMcTrack.h"
 #include "StMuDSTMaker/COMMON/StMuMcVertex.h"
+#include "StMuDSTMaker/COMMON/StMuFcsCollection.h"
+#include "StMuDSTMaker/COMMON/StMuFcsHit.h"
+#include "StMuDSTMaker/COMMON/StMuFcsCluster.h"
 
 #include "StTriggerUtilities/StTriggerSimuMaker.h"
 #include "StTriggerUtilities/Bemc/StBemcTriggerSimu.h"
@@ -70,6 +77,7 @@
 #include "StEmcRawMaker/StBemcTables.h"
 
 #include "StFmsDbMaker/StFmsDbMaker.h"
+#include "StFcsDbMaker/StFcsDb.h"
 
 #include "tables/St_mtdModuleToQTmap_Table.h"
 #include "tables/St_mtdQTSlewingCorr_Table.h"
@@ -94,6 +102,10 @@
 #include "StPicoEvent/StPicoBEmcSmdPHit.h"
 #include "StPicoEvent/StPicoETofHit.h"
 #include "StPicoEvent/StPicoETofPidTraits.h"
+#include "StPicoEvent/StPicoFwdTrack.h"
+#include "StPicoEvent/StPicoFwdVertex.h"
+#include "StPicoEvent/StPicoFcsHit.h"
+#include "StPicoEvent/StPicoFcsCluster.h"
 #include "StPicoEvent/StPicoMcVertex.h"
 #include "StPicoEvent/StPicoMcTrack.h"
 #include "StPicoEvent/StPicoArrays.h"
@@ -245,6 +257,10 @@ void StPicoDstMaker::streamerOff() {
   StPicoTrackCovMatrix::Class()->IgnoreTObjectStreamer();
   StPicoETofHit::Class()->IgnoreTObjectStreamer();
   StPicoETofPidTraits::Class()->IgnoreTObjectStreamer();
+  StPicoFwdTrack::Class()->IgnoreTObjectStreamer();
+  StPicoFwdVertex::Class()->IgnoreTObjectStreamer();
+  StPicoFcsHit::Class()->IgnoreTObjectStreamer();
+  StPicoFcsCluster::Class()->IgnoreTObjectStreamer();
   StPicoMcVertex::Class()->IgnoreTObjectStreamer();
   StPicoMcTrack::Class()->IgnoreTObjectStreamer();
 }
@@ -390,6 +406,11 @@ Int_t StPicoDstMaker::setVtxModeAttr(){
     return kStOK;
   }
 #endif /* !__TFG__VERSION__ */
+  else if (strcasecmp(SAttr("PicoVtxMode"), "PicoVtxless") == 0) {
+      setVtxMode( PicoVtxMode::Vtxless );
+      LOG_INFO << " PicoVtxless is being used " << endm;
+      return kStOK;
+    }
 
   return kStErr;
 }
@@ -893,6 +914,9 @@ Int_t StPicoDstMaker::MakeWrite() {
   fillEpdHits();
   fillBbcHits();
   fillETofHits();
+  fillFcsHits();
+  fillFcsClusters();
+  fillFwdTracks();
 
   // Could be a good idea to move this call to Init() or InitRun()
   StFmsDbMaker* fmsDbMaker = static_cast<StFmsDbMaker*>(GetMaker("fmsDb"));
@@ -904,7 +928,7 @@ Int_t StPicoDstMaker::MakeWrite() {
   mFmsFiller.fill(*mMuDst);
 
   if (Debug()) mPicoDst->printTracks();
-
+  
   mTTree->Fill();
   if ( isFromDaq ) {
 //    delete mEmcCollection;
@@ -1013,7 +1037,10 @@ void StPicoDstMaker::fillTracks() {
   // don't use StMuTrack::primary(), it returns primary tracks associated with
   // all vertices
   std::unordered_map<unsigned int, unsigned int> index2Primary;
-
+
+  if (mMuDst->primaryTracks() == nullptr) {
+    return;
+  }
   // Retrieve number of primary tracks
   Int_t nPrimarys = mMuDst->numberOfPrimaryTracks();
 
@@ -1064,7 +1091,7 @@ void StPicoDstMaker::fillTracks() {
     }
 
 #if defined (__TFG__VERSION__)
-    if ( StMuDst::dca3Dmax()>0 ) {
+    if ( StMuDst::dca3Dmax()>0) {
       // Cut large Dca
       THelixTrack t = dcaG->thelix();
       StThreeVectorD V(mMuDst->primaryVertex()->position());
@@ -1191,7 +1218,7 @@ void StPicoDstMaker::fillTracks() {
 
     // Calculate global momentum and position at point of DCA to the pVtx
     // at it is done in MuDst
-    if( dcaG ) {
+    if( dcaG && mMuDst->primaryVertex() ) {
       const StThreeVectorF &pvert =  mMuDst->primaryVertex()->position();
       Double_t vtx[3] = {pvert[0],pvert[1],pvert[2]};
       THelixTrack thelix =  dcaG->thelix();
@@ -1203,7 +1230,7 @@ void StPicoDstMaker::fillTracks() {
       picoTrk->setGlobalMomentum( momAtDca[0], momAtDca[1], momAtDca[2] );
       picoTrk->setOrigin( pos[0], pos[1], pos[2] );
     }
-    else {
+    else if ( mMuDst->primaryVertex() ){
       StPhysicalHelixD gHelix = dcaG->helix();
       gHelix.moveOrigin( gHelix.pathLength( mMuDst->primaryVertex()->position() ) );
       StThreeVectorF globMom =  gHelix.momentum( mBField * kilogauss );
@@ -1781,14 +1808,15 @@ void StPicoDstMaker::fillEvent() {
   picoEvent->setFillId( ev->runInfo().beamFillNumber( blue ) );
   picoEvent->setBField( ev->magneticField() );
   picoEvent->setTime( ev->eventInfo().time() );
-
   // Set primary vertex information
-  picoEvent->setPrimaryVertexPosition( ev->primaryVertexPosition().x(),
-				       ev->primaryVertexPosition().y(),
-				       ev->primaryVertexPosition().z() );
-  picoEvent->setPrimaryVertexPositionError( ev->primaryVertexErrors().x(),
-					    ev->primaryVertexErrors().y(),
-					    ev->primaryVertexErrors().z() );
+  if ( mMuDst->numberOfPrimaryVertices() > 0 ) {
+    picoEvent->setPrimaryVertexPosition( ev->primaryVertexPosition().x(),
+                ev->primaryVertexPosition().y(),
+                ev->primaryVertexPosition().z() );
+    picoEvent->setPrimaryVertexPositionError( ev->primaryVertexErrors().x(),
+                ev->primaryVertexErrors().y(),
+                ev->primaryVertexErrors().z() );
+  }
 
   if( StMuPrimaryVertex *pv = mMuDst->primaryVertex() ) {
     picoEvent->setPrimaryVertexRanking( pv->ranking() );
@@ -1828,7 +1856,11 @@ void StPicoDstMaker::fillEvent() {
 
   picoEvent->setGRefMult( ev->grefmult() );
   picoEvent->setNumberOfGlobalTracks( mMuDst->numberOfGlobalTracks() );
-  picoEvent->setNumberOfPrimaryTracks( mMuDst->numberOfPrimaryTracks() );
+  if (mMuDst->numberOfPrimaryVertices() > 0 && mMuDst->primaryTracks() ){ //calling mMudst->numberOfPrimaryTracks() will crash if there are no primary vertices
+    picoEvent->setNumberOfPrimaryTracks( mMuDst->numberOfPrimaryTracks() );
+  } else {
+    picoEvent->setNumberOfPrimaryTracks( 0 );
+  }
   picoEvent->setbTofTrayMultiplicity( ev->btofTrayMultiplicity() );
   picoEvent->setETofHitMultiplicity( ev->etofHitMultiplicity() );
   StMuETofHeader *header = mMuDst->etofHeader();
@@ -1927,7 +1959,7 @@ void StPicoDstMaker::fillEvent() {
   } //for(UInt_t iPMT=0; iPMT<bbc.numberOfPMTs(); ++iPMT)
 
     // Set bunch crossing ID
-  picoEvent->setBunchId( ev->l0Trigger().bunchCrossingId() );
+  picoEvent->setBunchId( ev->l0Trigger().bunchCrossingId7bit(ev->runNumber()) );
 }
 
 //_________________
@@ -2502,6 +2534,173 @@ void StPicoDstMaker::fillMtdHits() {
   } //while (triggerPos.size() > 0)
 }
 
+//_________________
+void StPicoDstMaker::fillFwdTracks() {
+
+  StEvent *evt = (StEvent *)GetDataSet("StEvent");
+  if ( evt ){
+    StFwdTrackCollection * evc = evt->fwdTrackCollection();
+    if ( !evc ){
+      LOG_ERROR << "null FwdTrackCollection" << endm;
+      return;
+    }
+
+    // fill a map of picodst FCS clusters between index and the cluster pointer
+    map<UShort_t, StPicoFcsCluster*> fcsClusterMap;
+    for ( auto iClu = 0; iClu < mPicoArrays[StPicoArrays::FcsCluster]->GetEntries(); iClu++ ){
+      StPicoFcsCluster * picoFcsCluster = (StPicoFcsCluster*)mPicoArrays[StPicoArrays::FcsCluster]->At(iClu);
+      fcsClusterMap[picoFcsCluster->index()] = picoFcsCluster; 
+    }
+
+    const StSPtrVecFwdTrack& evTracks = evc->tracks();
+    LOG_DEBUG << "Adding " << evc->numberOfTracks() << " StFwdTracks from StEvent to PicoDst" << endm; 
+    for ( size_t i = 0; i < evc->numberOfTracks(); i++ ){
+      StFwdTrack * evTrack = evTracks[i];
+      StPicoFwdTrack picoFwdTrack;
+      // Set the PicoDst attributes
+      picoFwdTrack.setMomentum( evTrack->momentum().x(), evTrack->momentum().y(), evTrack->momentum().z() );
+      picoFwdTrack.setNumberOfFitPoints( evTrack->numberOfFitPoints() * evTrack->charge() );
+      picoFwdTrack.setNumberOfSeedPoints( evTrack->numberOfSeedPoints() );
+      picoFwdTrack.setChi2( evTrack->chi2() );
+      picoFwdTrack.setPVal( evTrack->pval() );
+      picoFwdTrack.setDca( evTrack->dca().x(), evTrack->dca().y(), evTrack->dca().z() );
+      if ( evTrack->didFitConvergeFully())
+        picoFwdTrack.setStatus( 2 );
+      else if ( evTrack->didFitConverge())
+        picoFwdTrack.setStatus( 1 );
+      else
+        picoFwdTrack.setStatus( 0 );
+
+      picoFwdTrack.setMcTruth( evTrack->idTruth(), evTrack->qaTruth() );
+      picoFwdTrack.setVtxIndexRaw( evTrack->vertexIndexRaw() );
+      picoFwdTrack.setGlobalTrackIndex( evTrack->globalTrackIndex() );
+
+      // fill matched ecal and hcal clusters for the track
+      // ecal
+      for ( auto & cluster : evTrack->ecalClusters() ){
+        if ( mMapFcsIdPairIndex.count( make_pair( cluster->detectorId(), cluster->id() ) ) == 0 ){
+          LOG_WARN << "No FCS cluster found for ecal cluster with id " << cluster->id() << " and detectorId " << cluster->detectorId() << endm;
+          continue;
+        }
+        int index = mMapFcsIdPairIndex[ make_pair( cluster->detectorId(), cluster->id() ) ];
+        picoFwdTrack.addEcalCluster( index );
+      }
+      // hcal
+      for ( auto & cluster : evTrack->hcalClusters() ){
+        if ( mMapFcsIdPairIndex.count( make_pair( cluster->detectorId(), cluster->id() ) ) == 0 ){
+          LOG_WARN << "No FCS cluster found for hcal cluster with id " << cluster->id() << " and detectorId " << cluster->detectorId() << endm;
+          continue;
+        }
+        int index = mMapFcsIdPairIndex[ make_pair( cluster->detectorId(), cluster->id() ) ];
+        picoFwdTrack.addHcalCluster( index );
+      }
+
+      // Now set the projections for ECal and HCal
+      StFwdTrackProjection ecalProj = evTrack->getProjectionFor(kFcsWcalId);
+      StFwdTrackProjection hcalProj = evTrack->getProjectionFor(kFcsHcalId);
+      if ( !std::isnan( ecalProj.mXYZ.x() ) && !std::isnan( ecalProj.mXYZ.y() ) && !std::isnan( ecalProj.mXYZ.z() ) )
+        picoFwdTrack.setECalProjection( ecalProj.mXYZ.x(), ecalProj.mXYZ.y(), ecalProj.mXYZ.z() );
+      else 
+        picoFwdTrack.setECalProjection( -999., -999., -999. );
+      if ( !std::isnan( hcalProj.mXYZ.x() ) && !std::isnan( hcalProj.mXYZ.y() ) && !std::isnan( hcalProj.mXYZ.z() ) )
+        picoFwdTrack.setHCalProjection( hcalProj.mXYZ.x(), hcalProj.mXYZ.y(), hcalProj.mXYZ.z() );
+      else
+        picoFwdTrack.setHCalProjection( -999., -999., -999. );
+
+      int counter = mPicoArrays[StPicoArrays::FwdTrack]->GetEntries();
+      picoFwdTrack.setId( counter );
+      new((*(mPicoArrays[StPicoArrays::FwdTrack]))[counter]) StPicoFwdTrack(picoFwdTrack);
+    }
+  } else {
+    LOG_DEBUG << "Cannot get Fwd Tracks from StEvent" << endm;
+  }
+
+  // Fill FwdVertex also
+  // get primary vertex from StEvent
+  LOG_DEBUG << "There are " << evt->numberOfPrimaryVertices() << " primary vertices from StEvent" << endm; 
+  for ( size_t i = 0; i < evt->numberOfPrimaryVertices(); i++ ){
+    StPrimaryVertex * evVertex = evt->primaryVertex(i);
+    if (!evVertex || !evVertex->isFwdVtx()) continue; // only save valid FwdVertex
+    StPicoFwdVertex picoFwdVertex;
+    // Set the PicoDst attributes
+    picoFwdVertex.setPosition( evVertex->position().x(), evVertex->position().y(), evVertex->position().z() );
+    LOG_DEBUG << "Adding StPicoFwdVertex from StPrimaryVertex" << endm;
+    picoFwdVertex.setChi2( evVertex->chiSquared() );
+    picoFwdVertex.setNumberOfTracks( evVertex->numTracksUsedInFinder() );
+    int counter = mPicoArrays[StPicoArrays::FwdVertex]->GetEntries();
+    new((*(mPicoArrays[StPicoArrays::FwdVertex]))[counter]) StPicoFwdVertex(picoFwdVertex);
+  } // for each primary vertex
+
+} //fillFwdTracks
+
+//_________________
+void StPicoDstMaker::fillFcsClusters() {
+  StFcsDb* fcsDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
+  if( !fcsDb ) {
+      LOG_ERROR << "Cannot get StFcsDb object" << endm;
+      return;
+  }
+  StMuFcsCollection * muFcs = mMuDst->muFcsCollection();
+  if ( !muFcs ) {
+    LOG_ERROR << "Cannot get Fcs Collection from MuDst" << endm;
+    return;
+  }
+  TIter next(muFcs->getClusterArray());
+  StMuFcsCluster* muCluster(NULL);
+  while ((muCluster = static_cast<StMuFcsCluster*>(next()))) {
+    int counter = mPicoArrays[StPicoArrays::FcsCluster]->GetEntries();
+    StPicoFcsCluster picoFcsCluster;
+    picoFcsCluster.setId(muCluster->id());
+    picoFcsCluster.setIndex(counter);
+    picoFcsCluster.setDetectorId(muCluster->detectorId());
+    picoFcsCluster.setCategory(muCluster->category());
+    picoFcsCluster.setNTowers(muCluster->nTowers());
+    picoFcsCluster.setEnergy(muCluster->energy());
+    picoFcsCluster.setX(muCluster->x());
+    picoFcsCluster.setY(muCluster->y());
+    picoFcsCluster.setSigmaMin(muCluster->sigmaMin());
+    picoFcsCluster.setSigmaMax(muCluster->sigmaMax());
+    picoFcsCluster.setTheta(muCluster->theta());
+    picoFcsCluster.setChi2Ndf1Photon(muCluster->chi2Ndf1Photon());
+    picoFcsCluster.setChi2Ndf2Photon(muCluster->chi2Ndf2Photon());
+    double zVertex=picoDst()->event()->primaryVertex().z();
+    StThreeVectorD xyz=fcsDb->getStarXYZfromColumnRow(muCluster->detectorId(), muCluster->x(), muCluster->y());    
+    StLorentzVectorD lv = fcsDb->getLorentzVector(xyz, muCluster->energy(), zVertex);
+    picoFcsCluster.setFourMomentum(lv.px(),lv.py(),lv.pz(),lv.e());
+    new((*(mPicoArrays[StPicoArrays::FcsCluster]))[counter]) StPicoFcsCluster(picoFcsCluster);
+    mMapFcsIdPairIndex[ make_pair( muCluster->detectorId(), muCluster->id() ) ] = counter;
+  }
+  LOG_INFO << "StPicoDstMaker::fillFcsClusters filled " << mPicoArrays[StPicoArrays::FcsCluster]->GetEntries() << endm;
+}//fillFcsClusters
+
+//_________________
+void StPicoDstMaker::fillFcsHits() {
+  StFcsDb* fcsDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
+  if( !fcsDb ) {
+      LOG_ERROR << "Cannot get StFcsDb object" << endm;
+      return;
+  }
+  StMuFcsCollection * muFcs = mMuDst->muFcsCollection();
+  if ( !muFcs ) {
+    LOG_ERROR << "Cannot get Fcs Collection from MuDst" << endm;
+    return;
+  }
+  TIter next(muFcs->getHitArray());
+  StMuFcsHit* muHit(NULL);
+  while ((muHit = static_cast<StMuFcsHit*>(next()))) {
+    int counter = mPicoArrays[StPicoArrays::FcsHit]->GetEntries();
+    StPicoFcsHit picoFcsHit;
+    picoFcsHit.setDetectorId(muHit->detectorId());
+    picoFcsHit.setId(muHit->id());
+    double zVertex=picoDst()->event()->primaryVertex().z();
+    StThreeVectorD xyz = fcsDb->getStarXYZ(muHit->detectorId(), muHit->id());
+    StLorentzVectorD lv = fcsDb->getLorentzVector(xyz, muHit->energy(), zVertex);
+    picoFcsHit.setFourMomentum(lv.px(),lv.py(),lv.pz(),lv.e());
+    new((*(mPicoArrays[StPicoArrays::FcsHit]))[counter]) StPicoFcsHit(picoFcsHit);
+  }
+  LOG_INFO << "StPicoDstMaker::fillFcsHits filled " << mPicoArrays[StPicoArrays::FcsHit]->GetEntries() << endm;
+}//fillFcsHit
+
 #if !defined (__TFG__VERSION__)
 
 /**
@@ -2523,6 +2722,9 @@ bool StPicoDstMaker::selectVertex() {
   StMuPrimaryVertex* selectedVertex = nullptr;
 
   // Switch between modes: Default and Vpd (VpdOrDefault)
+  if ( mVtxMode == PicoVtxMode::Vtxless ) {
+    return true;
+  }
   // Default takes the first primary vertex, meanwhile
   // Vpd
   if ( mVtxMode == PicoVtxMode::Default ) {

StRoot/StPicoDstMaker/StPicoDstMaker.h

@@ -90,7 +90,8 @@ class StPicoDstMaker : public StMaker {
   /// \par 3 VpdOrDefault
   /// \par 4 Mtd
   /// \par 5 FXT
-  enum PicoVtxMode {NotSet=0, Default=1, Vpd=2, VpdOrDefault=3, Mtd=4, FXT=5};
+  /// \par 6 VTXLess
+  enum PicoVtxMode {NotSet=0, Default=1, Vpd=2, VpdOrDefault=3, Mtd=4, FXT=5, Vtxless=6};
 #endif /* ! __TFG__VERSION__ */
   /// Write or not write covariance matrix: 0-skip, 1-write
   enum PicoCovMtxMode {Skip=0, Write=1};
@@ -225,6 +226,12 @@ class StPicoDstMaker : public StMaker {
   void fillBbcHits();
   /// Fill ETOF information
   void fillETofHits();
+  /// Fill Fwd Track information
+  void fillFwdTracks();
+  /// Fill FcsHits information
+  void fillFcsHits();
+  /// Fill FcsClusters information
+  void fillFcsClusters();
   /// Fill MC vertex information
   void fillMcVertices();
   /// Fill MC track information
@@ -287,6 +294,9 @@ class StPicoDstMaker : public StMaker {
   /// Pointer to the array of BEMC tower hits
   StEmcRawHit*     mEmcIndex[4800];
 
+  /// Map of StPicoFcsCluster detId + Id to index in the array for track matching
+  std::map<std::pair<int, int>, int> mMapFcsIdPairIndex;
+
   /// Magnetic field of the current event
   Float_t    mBField;