Nikita Ermakov 1 year ago
parent
commit
e79269d7fc
100 changed files with 30489 additions and 2990 deletions
  1. 462 275
      StRoot/StFemtoDstMaker/StFemtoDstMaker.cxx
  2. 22 1
      StRoot/StFemtoDstMaker/StFemtoDstMaker.h
  3. 9 1
      StRoot/StFemtoDstMaker/StFemtoDstQAMaker.cxx
  4. 3 0
      StRoot/StFemtoDstMaker/StFemtoDstQAMaker.h
  5. 31 26
      StRoot/StFemtoDstMaker/StFemtoEvent.h
  6. 2 0
      StRoot/StFemtoDstMaker/StHbtFemtoDstReader.cxx
  7. 40 20
      StRoot/StFemtoDstMaker/StHbtO97DstReader.cxx
  8. 11 7
      StRoot/StFemtoDstMaker/StHbtO97DstReader.h
  9. 34 21
      StRoot/StFemtoDstMaker/StO97DstMaker.cxx
  10. 1 0
      StRoot/StFemtoDstMaker/StO97DstMaker.h
  11. 305 252
      StRoot/StFemtoDstMaker/StO97DstQAMaker.cxx
  12. 31 12
      StRoot/StFemtoDstMaker/StO97DstQAMaker.h
  13. 3 2
      StRoot/StFemtoDstMaker/StO97Event.cxx
  14. 0 1
      StRoot/StFemtoDstMaker/StO97Event.h
  15. 21 26
      StRoot/StFemtoDstMaker/StO97Track.cxx
  16. 17 3
      StRoot/StFemtoDstMaker/StO97Track.h
  17. 3797 0
      StRoot/StFemtoDstMaker/pdg_table.txt
  18. 592 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstInclusiveSelector.cxx
  19. 228 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstInclusiveSelector.h
  20. 689 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker.cxx
  21. 218 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker.h
  22. 669 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker_HFT.cxx
  23. 206 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker_HFT.h
  24. 535 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker_StarGen.cxx
  25. 115 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstMaker_StarGen.h
  26. 508 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPicoRun11Maker.cxx
  27. 207 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPicoRun11Maker.h
  28. 510 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPicoRun12Maker.cxx
  29. 207 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPicoRun12Maker.h
  30. 470 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPythia6Maker.cxx
  31. 107 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstPythia6Maker.h
  32. 481 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstQAMaker.cxx
  33. 118 0
      StRoot/StFemtoDstMakerNORP/StFemtoDstQAMaker.h
  34. 193 0
      StRoot/StFemtoDstMakerNORP/StFemtoEvent.cxx
  35. 114 0
      StRoot/StFemtoDstMakerNORP/StFemtoEvent.h
  36. 337 0
      StRoot/StFemtoDstMakerNORP/StFemtoTrack.cxx
  37. 133 0
      StRoot/StFemtoDstMakerNORP/StFemtoTrack.h
  38. 393 0
      StRoot/StFemtoDstMakerNORP/StHbtFemtoDstReader.cxx
  39. 76 0
      StRoot/StFemtoDstMakerNORP/StHbtFemtoDstReader.h
  40. 406 0
      StRoot/StFemtoDstMakerNORP/StHbtO97DstReader.cxx
  41. 68 0
      StRoot/StFemtoDstMakerNORP/StHbtO97DstReader.h
  42. 286 0
      StRoot/StFemtoDstMakerNORP/StO97DstMaker.cxx
  43. 82 0
      StRoot/StFemtoDstMakerNORP/StO97DstMaker.h
  44. 359 0
      StRoot/StFemtoDstMakerNORP/StO97DstQAMaker.cxx
  45. 105 0
      StRoot/StFemtoDstMakerNORP/StO97DstQAMaker.h
  46. 90 0
      StRoot/StFemtoDstMakerNORP/StO97Event.cxx
  47. 63 0
      StRoot/StFemtoDstMakerNORP/StO97Event.h
  48. 451 0
      StRoot/StFemtoDstMakerNORP/StO97FemtoDstMaker.cxx
  49. 129 0
      StRoot/StFemtoDstMakerNORP/StO97FemtoDstMaker.h
  50. 142 0
      StRoot/StFemtoDstMakerNORP/StO97Track.cxx
  51. 99 0
      StRoot/StFemtoDstMakerNORP/StO97Track.h
  52. 3797 0
      StRoot/StFemtoDstMakerNORP/pdg_table.txt
  53. 20 16
      StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt.cxx
  54. 192 0
      StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_th.cxx
  55. 67 0
      StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_th.h
  56. 144 148
      StRoot/StHbtMaker/Cut/gregsEventCut.cxx
  57. 9 1
      StRoot/StHbtMaker/Cut/gregsEventCut.h
  58. 2 1
      StRoot/StHbtMaker/Cut/gregsPionTrackCut.cxx
  59. 286 0
      StRoot/StHbtMaker/Cut/gregsProtonTrackCut.cxx
  60. 103 0
      StRoot/StHbtMaker/Cut/gregsProtonTrackCut.h
  61. 5 2
      StRoot/StHbtMaker/Infrastructure/StHbtEvent.hh
  62. 10 0
      StRoot/StHbtMaker/Infrastructure/StHbtTrack.cc
  63. 9 1
      StRoot/StHbtMaker/Infrastructure/StHbtTrack.hh
  64. 438 452
      StRoot/StHbtStarGenReader/StHbtStarGenReader.cxx
  65. 1 4
      StRoot/StHbtStarGenReader/StHbtStarGenReader.h
  66. 339 0
      StRoot/StQGSM/StQGSMdstMaker.cxx
  67. 73 0
      StRoot/StQGSM/StQGSMdstMaker.h
  68. 356 0
      StRoot/StQGSM/StQGSMdstQAMaker.cxx
  69. 105 0
      StRoot/StQGSM/StQGSMdstQAMaker.h
  70. 2 0
      StRoot/StSphericityStarGen/StSphericity.cxx
  71. 2 2
      macros/HBT/expKaonHBT_auau200_3D.C
  72. 0 530
      macros/HBT/pp/expPionHBT_pp200_y2012_1D.C
  73. 15 5
      macros/HBT/pp/nomult/expKaonHBT_pp200_y2012_1D.C
  74. 46 39
      macros/HBT/pp/nomult/expPionHBT_pp200_y2012_1D.C
  75. 517 0
      macros/HBT/pp/nomult/expUnlikeKaonHBT_pp200_y2012_1D.C
  76. 517 0
      macros/HBT/pp/nomult/expUnlikePionHBT_pp200_y2012_1D.C
  77. 476 0
      macros/HBT/pp/nomult/sim97KaonHBT_pp200n510_1D.C
  78. 478 0
      macros/HBT/pp/nomult/sim97PionHBT_pp200n510_1D.C
  79. 479 0
      macros/HBT/pp/nomult/sim97ProtonHBT_pp200n510_1D.C
  80. 59 51
      macros/HBT/pp/nomult/simPionHBT_pp200n510_1D.C
  81. 506 0
      macros/HBT/pp/nomult/simProtonHBT_pp200n510_1D.C
  82. 506 0
      macros/HBT/pp/nomult/simSGPionHBT_pp200n510_1D.C
  83. 494 0
      macros/HBT/pp/nomult/simUnlikePionHBT_pp200_y2012_1D.C
  84. 179 178
      macros/HBT/pp/expKaonHBT_pp200_y2012_1D.C
  85. 172 172
      macros/HBT/pp/expKaonHBT_pp510_y2011_1D.C
  86. 534 0
      macros/HBT/pp/nomult/tof/expPionHBT_pp200_y2012_1D_TOF.C
  87. 178 177
      macros/HBT/pp/expPionHBT_pp510_y2011_1D.C
  88. 509 0
      macros/HBT/pp/nomult/tof/simKaonHBT_pp200n510_1D_TOF.C
  89. 516 0
      macros/HBT/pp/nomult/tof/simPionHBT_pp200n510_1D_TOF.C
  90. 531 0
      macros/HBT/pp/nomult/tpc/expKaonHBT_pp200_y2012_1D_TPC.C
  91. 531 0
      macros/HBT/pp/nomult/tpc/expKaonHBT_pp510_y2011_1D_TPC.C
  92. 541 0
      macros/HBT/pp/nomult/tpc/expPionHBT_pp200_y2012_1D_TPC.C
  93. 533 0
      macros/HBT/pp/nomult/tpc/expPionHBT_pp510_y2011_1D_TPC.C
  94. 509 0
      macros/HBT/pp/nomult/tpc/simKaonHBT_pp200n510_1D_TPC.C
  95. 517 0
      macros/HBT/pp/nomult/tpc/simPionHBT_pp200n510_1D_TPC.C
  96. 162 191
      macros/HBT/pp/unlike/expKaonHBT_pp200_y2012_1D_unlike.C
  97. 153 184
      macros/HBT/pp/unlike/expKaonHBT_pp510_y2011_1D_unlike.C
  98. 534 0
      macros/HBT/pp/nomult/tpc/tof/expPionHBT_pp200_y2012_1D_TOF.C
  99. 162 189
      macros/HBT/pp/unlike/expPionHBT_pp510_y2011_1D_unlike.C
  100. 0 0
      macros/HBT/pp/nomult/tpc/tof/simKaonHBT_pp200n510_1D_TOF.C

+ 462 - 275
StRoot/StFemtoDstMaker/StFemtoDstMaker.cxx

@@ -13,6 +13,7 @@ ClassImp(StFemtoDstMaker)
 
 //_________________
 StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
+                                 const char *calibPath,                                 
                                  const Char_t *oFileName) {
 
   std::cout << "StFemtoDstMaker::StFemtoDstMaker - Creating an instance...";
@@ -43,6 +44,8 @@ StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
   mIsGoodTrack = false;
   mCurrentTrigger = 0;
 
+  mCalibPath = calibPath;
+
   //
   // Default cut values
   //
@@ -96,7 +99,8 @@ StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
 
 //_________________
 StFemtoDstMaker::~StFemtoDstMaker() {
-    delete matrix;
+  delete mQVMaker;
+  delete matrix;
 }
 
 //_________________
@@ -122,11 +126,14 @@ Int_t StFemtoDstMaker::Init() {
       mRefMultCorrUtil = new StRefMultCorr("refmult");
       std::cout << "\t[DONE]" << endl;
     }
+
+    // initialize reaction plane maker
+    mQVMaker = new QVMaker(4 /* do not fill hEP histograms just calculate Psi */);
   }
 
   mNBytes = 0;
   std::cout << "StFemtoDstMaker::Init - Initialization has been finished"
-      << std::endl;
+            << std::endl;
   return StMaker::Init();
 }
 
@@ -150,7 +157,7 @@ Int_t StFemtoDstMaker::Make() {
     mMuDstMaker = (StMuDstMaker*)GetMaker("MuDst");
     if (!mMuDstMaker) {
       LOG_ERROR << "StFemtoDstMaker::Make [ERROR] - Cannot find StMuDstMaker"
-          << std::endl;
+                << std::endl;
       return kStOk;
     }
   }
@@ -160,12 +167,12 @@ Int_t StFemtoDstMaker::Make() {
   //
   mMuDst = mMuDstMaker->muDst();
   if (!mMuDst) {
-      mMuDst = (StMuDst*)GetInputDS("MuDst");
-      if (!mMuDst) {
-          gMessMgr->Warning() << "StFemtoDstMaker::Make [WARNING] - No MuDst has been found" 
-              << endm;
-          return kStOk;
-      }
+    mMuDst = (StMuDst*)GetInputDS("MuDst");
+    if (!mMuDst) {
+      gMessMgr->Warning() << "StFemtoDstMaker::Make [WARNING] - No MuDst has been found" 
+                          << endm;
+      return kStOk;
+    }
   }
 
   //
@@ -174,8 +181,8 @@ Int_t StFemtoDstMaker::Make() {
   mMuEvent = (StMuEvent*)mMuDst->event();
 
   if(mMuEvent->refMult() < 0) { // multiplicity cannot be negative
-      mRefMultFail++;
-      return kStOk;
+    mRefMultFail++;
+    return kStOk;
   }
 
   Int_t mNVertices = mMuDst->numberOfPrimaryVertices();
@@ -192,8 +199,8 @@ Int_t StFemtoDstMaker::Make() {
 
   if (mPrimVertex->ranking() <= 0) // positive ranking only
   {
-      mRankingFail++;
-      return kStOk;
+    mRankingFail++;
+    return kStOk;
   }
 
   //
@@ -203,14 +210,14 @@ Int_t StFemtoDstMaker::Make() {
      mPrimVertex->position().y() == 0 &&
      mPrimVertex->position().z() == 0)
   {
-      mVtxNullFail++;
-      return kStOk;
+    mVtxNullFail++;
+    return kStOk;
   }
 
   if(mNPrimTracks < 0 || mNPrimTracks > 10000) // reasonable amount of tracks
   {
-      mNPrimVtxFail++;
-      return kStOk;
+    mNPrimVtxFail++;
+    return kStOk;
   }
 
   //
@@ -225,31 +232,33 @@ Int_t StFemtoDstMaker::Make() {
   Double_t mMassSqr;
   Int_t   mPrimVertIndex = -999;
   Float_t mRanking = -999.;
+  Bool_t mGoodRPevent = false;
+  
 
   if(mCollisionType == false) { // pp collisions means false state
-      mCent16 = -1;
-      mCentWeight = -999.;
+    mCent16 = -1;
+    mCentWeight = -999.;
   }
 
   //Calculate the amount of the pile-up vertices
   int mNVtxInBunchCross = 0;
   for (unsigned int i = 0; i < mMuDst->numberOfPrimaryVertices(); i++)
   {
-      StMuPrimaryVertex *primVtx = mMuDst->primaryVertex(i);
-      if(!primVtx || primVtx->ranking() <= 0.)   continue;
+    StMuPrimaryVertex *primVtx = mMuDst->primaryVertex(i);
+    if(!primVtx || primVtx->ranking() <= 0.)   continue;
 
-      if(TMath::Abs(primVtx->position().x()) < 1e-5 &&
-         TMath::Abs(primVtx->position().y()) < 1e-5 &&
-         TMath::Abs(primVtx->position().z()) < 1e-5)     continue;
+    if(TMath::Abs(primVtx->position().x()) < 1e-5 &&
+       TMath::Abs(primVtx->position().y()) < 1e-5 &&
+       TMath::Abs(primVtx->position().z()) < 1e-5)     continue;
 
-      if(TMath::Abs( primVtx->position().z() - mTofHeader->vpdVz()) < 5. )   mNVtxInBunchCross++;
+    if(TMath::Abs( primVtx->position().z() - mTofHeader->vpdVz()) < 5. )   mNVtxInBunchCross++;
   }
 
   //There should be only one primary vertex that caused the trigger
   if(mNVtxInBunchCross > 1)
   {
-      mAntiPileupFail++;
-      return kStOk;
+    mAntiPileupFail++;
+    return kStOk;
   }
 
   mEventIsGood = false;
@@ -261,14 +270,14 @@ Int_t StFemtoDstMaker::Make() {
 
   if(!AcceptTrigger(mMuEvent)) // if trigger is not found
   {
-      mAcceptTriggerFail++;
-      return kStOk;
+    mAcceptTriggerFail++;
+    return kStOk;
   }
 
   if(!AcceptPrimVtx(mVertPosition, mVpdVz))
   {
-      mAcceptPrimVtxFail++;
-      return kStOk;
+    mAcceptPrimVtxFail++;
+    return kStOk;
   }
 
   mEventIsGood = true;
@@ -277,21 +286,34 @@ Int_t StFemtoDstMaker::Make() {
   //
   if(mCollisionType == true) // AuAu collisions means true state
   {
+    int mass1 = mMuEvent->runInfo().beamMassNumber(StBeamDirection::east);
+    int mass2 = mMuEvent->runInfo().beamMassNumber(StBeamDirection::west);
+    
+    // if Au+Au
+    if (mass1 == 197 && mass2 == 197) {
+
       mRefMultCorrUtil->init(mMuEvent->runId());
       if(mAuAuCorrectZdc == true) {
-          mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
-                                      mVertPosition.z(),
-                                      mMuEvent->runInfo().zdcCoincidenceRate());
+        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
+                                    mVertPosition.z(),
+                                    mMuEvent->runInfo().zdcCoincidenceRate());
       }
       else {
-          mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
-                                      mVertPosition.z());
+        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
+                                    mVertPosition.z());
       }
       mCent16 = mRefMultCorrUtil->getCentralityBin16();
       if(mCent16 < 0) // some AuAu collisions have -1 value
-          return kStOk;
+        return kStOk;
       mCentWeight = mRefMultCorrUtil->getWeight();
       mRefMultCorrVal = mRefMultCorrUtil->getRefMultCorr();
+    }
+    else { // Assume that it is Cu+Au case
+      mCent16 = GetCentralityBin(mMuEvent->refMult(), mVertPosition.z());
+    }
+
+    if (mCalibPath && LoadEPCalibParam(mMuEvent->runId()))
+      mGoodRPevent = true;
   } //if(mCollisionType == true)
 
   mNEventsPassed++;
@@ -322,7 +344,7 @@ Int_t StFemtoDstMaker::Make() {
   StFemtoTrack *mFTrack = NULL;
 
   //
-  // Calculate sphericity
+  // Calculate transverse sphericity
   //
   float sph          = -999.;
   float pt_full      = 0.0;
@@ -330,23 +352,23 @@ Int_t StFemtoDstMaker::Make() {
   matrix->Zero();
   for (int i = 0; i < mNPrimTracks; i++)
   {
-      mPrimTrack = mMuDst->primaryTracks(i);
-      if (mPrimTrack->vertexIndex() != 0)   return kStOk;
-
-      float pt = mPrimTrack->pt();
-      const StThreeVectorF &p = mPrimTrack->p();
-      if (fabs(mPrimTrack->eta()) < 1.0 && pt > mPtCut)
-      {
-          float px = p.x();
-          float py = p.y();
-
-          (*matrix)(0, 0) += px*px/pt;
-          (*matrix)(1, 1) += py*py/pt;
-          (*matrix)(0, 1) += px*py/pt;
-          (*matrix)(1, 0) += px*py/pt;
-
-          pt_full += pt;
-      }
+    mPrimTrack = mMuDst->primaryTracks(i);
+    if (mPrimTrack->vertexIndex() != 0)   return kStOk;
+
+    float pt = mPrimTrack->pt();
+    const StThreeVectorF &p = mPrimTrack->p();
+    if (fabs(mPrimTrack->eta()) < 1.0 && pt > mPtCut)
+    {
+      float px = p.x();
+      float py = p.y();
+
+      (*matrix)(0, 0) += px*px/pt;
+      (*matrix)(1, 1) += py*py/pt;
+      (*matrix)(0, 1) += px*py/pt;
+      (*matrix)(1, 0) += px*py/pt;
+
+      pt_full += pt;
+    }
   }
 
   *matrix *= 1./pt_full;
@@ -356,102 +378,123 @@ Int_t StFemtoDstMaker::Make() {
   mFemtoEvent->SetSphericity(sph);
 
   //
+  // Calculate reaction plane
+  //
+  float rpE = -999., rpW = -999.;
+  if (mGoodRPevent) {
+    int qZVtx = int((mPrimVertex->position().z() + 30.)/60. * NqvZvtx /* see ConstVar.h */);
+
+    if (mCent16 >= 0 && mCent16 < NqvCent &&
+        qZVtx >= 0 && qZVtx < NqvZvtx) {
+      
+      // SMD
+      QV Q_SMD[NsubSMD];
+      CalcQvSMD(mMuEvent, Q_SMD, mCent16, qZVtx);
+      rpE = Q_SMD[0].get_Psi();
+      rpW = Q_SMD[1].get_Psi();
+    }
+  }
+
+  mFemtoEvent->SetRPeast(rpE);
+  mFemtoEvent->SetRPwest(rpW);
+  
+  //
   // Loop over primary tracks
   //
   for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++)
   {
-      mPrimTrack = mMuDst->primaryTracks(iTrk);
-      int idxGlob = mPrimTrack->index2Global();
-      if (idxGlob < 0)   continue;
-      mGlobTrack = mMuDst->globalTracks(idxGlob);
-
-      if (AcceptRefMult2(mPrimTrack))
-      {
-          refMult2++;
-          if(mPrimTrack->charge() > 0)    refMult2Pos++;
-      }
+    mPrimTrack = mMuDst->primaryTracks(iTrk);
+    int idxGlob = mPrimTrack->index2Global();
+    if (idxGlob < 0)   continue;
+    mGlobTrack = mMuDst->globalTracks(idxGlob);
+
+    if (AcceptRefMult2(mPrimTrack))
+    {
+      refMult2++;
+      if(mPrimTrack->charge() > 0)    refMult2Pos++;
+    }
 
-      if (!AcceptTrack(mPrimTrack, iVert))   continue;
-
-      Bool_t mIsTofTrack = IsTofTrack(mGlobTrack);
-      if(mIsTofTrack) {
-          mBeta = mGlobTrack->btofPidTraits().beta();
-          mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
-          //
-          // Next lines are cheat. But we should remove garbage
-          //
-          if(mMassSqr < -0.02) continue;
-          if(mPrimTrack->momentum().mag() > 0.5) {
-              if(mMassSqr > 0.06 && mMassSqr < 0.18)   continue;
-              if(mMassSqr > 0.35 && mMassSqr < 0.8)    continue;
-              if(mMassSqr > 1.)                        continue;
-          }
-      } //if(mIsTofTrack)
-      else {
-          mBeta = -999.;
-          mMassSqr = -999.;
+    if (!AcceptTrack(mPrimTrack, iVert))   continue;
+
+    Bool_t mIsTofTrack = IsTofTrack(mGlobTrack);
+    if(mIsTofTrack) {
+      mBeta = mGlobTrack->btofPidTraits().beta();
+      mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
+      //
+      // Next lines are cheat. But we should remove garbage
+      //
+      if(mMassSqr < -0.02) continue;
+      if(mPrimTrack->momentum().mag() > 0.5) {
+        if(mMassSqr > 0.06 && mMassSqr < 0.18)   continue;
+        if(mMassSqr > 0.35 && mMassSqr < 0.8)    continue;
+        if(mMassSqr > 1.)                        continue;
       }
+    } //if(mIsTofTrack)
+    else {
+      mBeta = -999.;
+      mMassSqr = -999.;
+    }
 
-      Bool_t mIsTpcPion   = false;
-      Bool_t mIsTpcKaon   = false;
-      Bool_t mIsTpcProton = false;
-      Bool_t mIsTofPion   = false;
-      Bool_t mIsTofKaon   = false;
-      Bool_t mIsTofProton = false;
+    Bool_t mIsTpcPion   = false;
+    Bool_t mIsTpcKaon   = false;
+    Bool_t mIsTpcProton = false;
+    Bool_t mIsTofPion   = false;
+    Bool_t mIsTofKaon   = false;
+    Bool_t mIsTofProton = false;
 
-      if(mRemovePions) {
-          mIsTpcPion = IsTpcPion(mPrimTrack);
-          mIsTofPion = IsTofPion(mGlobTrack);
+    if(mRemovePions) {
+      mIsTpcPion = IsTpcPion(mPrimTrack);
+      mIsTofPion = IsTofPion(mGlobTrack);
 
-          if(mIsTpcPion==true || mIsTofPion==true)
-              continue;
-      }
+      if(mIsTpcPion==true || mIsTofPion==true)
+        continue;
+    }
 
-      if(mRemoveKaons) {
-          mIsTpcKaon = IsTpcKaon(mPrimTrack);
-          mIsTofKaon = IsTofKaon(mGlobTrack);
+    if(mRemoveKaons) {
+      mIsTpcKaon = IsTpcKaon(mPrimTrack);
+      mIsTofKaon = IsTofKaon(mGlobTrack);
 
-          if(mIsTpcKaon==true || mIsTofKaon==true)
-              continue;
-      }
+      if(mIsTpcKaon==true || mIsTofKaon==true)
+        continue;
+    }
 
-      if(mRemoveProtons) {
-          mIsTpcProton = IsTpcProton(mPrimTrack);
-          mIsTofProton = IsTofProton(mGlobTrack);
+    if(mRemoveProtons) {
+      mIsTpcProton = IsTpcProton(mPrimTrack);
+      mIsTofProton = IsTofProton(mGlobTrack);
 
-          if(mIsTpcProton==true || mIsTofProton==true)
-              continue;
-      }
+      if(mIsTpcProton==true || mIsTofProton==true)
+        continue;
+    }
 
-      mFTrack = mFemtoEvent->AddFemtoTrack();
-
-      mFTrack->SetId( mPrimTrack->id() );
-      mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
-      mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss() );
-      //mFTrack->SetNHitsDedx( mPrimTrack->nHitsDedx() );
-      mFTrack->SetNSigmaElectron( mPrimTrack->nSigmaElectron() );
-      mFTrack->SetNSigmaPion( mPrimTrack->nSigmaPion() );
-      mFTrack->SetNSigmaKaon( mPrimTrack->nSigmaKaon() );
-      mFTrack->SetNSigmaProton( mPrimTrack->nSigmaProton() );
-      mFTrack->SetDedx( mPrimTrack->dEdx() );
-      mFTrack->SetMap0( mGlobTrack->topologyMap().data(0) );
-      mFTrack->SetMap1( mGlobTrack->topologyMap().data(1) );
-      mFTrack->SetP( mPrimTrack->momentum().x(), 
-                    mPrimTrack->momentum().y(), 
-                    mPrimTrack->momentum().z() );
-      mFTrack->SetDCAxGlobal( mPrimTrack->dcaGlobal().x() );
-      mFTrack->SetDCAyGlobal( mPrimTrack->dcaGlobal().y() );
-      mFTrack->SetDCAzGlobal( mPrimTrack->dcaGlobal().z() );
-      mFTrack->SetGlobalP( mGlobTrack->p().x(),
-                          mGlobTrack->p().y(),
-                          mGlobTrack->p().z() );
-
-      if(mIsTofTrack) {
-          mFTrack->SetBeta(mBeta);
-      }
-      else {
-          mFTrack->SetBeta(-999.);
-      }
+    mFTrack = mFemtoEvent->AddFemtoTrack();
+
+    mFTrack->SetId( mPrimTrack->id() );
+    mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
+    mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss() );
+    //mFTrack->SetNHitsDedx( mPrimTrack->nHitsDedx() );
+    mFTrack->SetNSigmaElectron( mPrimTrack->nSigmaElectron() );
+    mFTrack->SetNSigmaPion( mPrimTrack->nSigmaPion() );
+    mFTrack->SetNSigmaKaon( mPrimTrack->nSigmaKaon() );
+    mFTrack->SetNSigmaProton( mPrimTrack->nSigmaProton() );
+    mFTrack->SetDedx( mPrimTrack->dEdx() );
+    mFTrack->SetMap0( mGlobTrack->topologyMap().data(0) );
+    mFTrack->SetMap1( mGlobTrack->topologyMap().data(1) );
+    mFTrack->SetP( mPrimTrack->momentum().x(), 
+                   mPrimTrack->momentum().y(), 
+                   mPrimTrack->momentum().z() );
+    mFTrack->SetDCAxGlobal( mPrimTrack->dcaGlobal().x() );
+    mFTrack->SetDCAyGlobal( mPrimTrack->dcaGlobal().y() );
+    mFTrack->SetDCAzGlobal( mPrimTrack->dcaGlobal().z() );
+    mFTrack->SetGlobalP( mGlobTrack->p().x(),
+                         mGlobTrack->p().y(),
+                         mGlobTrack->p().z() );
+
+    if(mIsTofTrack) {
+      mFTrack->SetBeta(mBeta);
+    }
+    else {
+      mFTrack->SetBeta(-999.);
+    }
 
   } //for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++)
 
@@ -459,8 +502,8 @@ Int_t StFemtoDstMaker::Make() {
   mFemtoEvent->SetRefMult2Pos(refMult2Pos);
 
   if(mEventIsGood) {
-      mNBytes += mTree->Fill();
-      mFemtoEvent->Clear();
+    mNBytes += mTree->Fill();
+    mFemtoEvent->Clear();
   }
 
 
@@ -469,221 +512,365 @@ Int_t StFemtoDstMaker::Make() {
 
 //_________________
 Bool_t StFemtoDstMaker::AcceptTrigger(StMuEvent *muEvent) {
-    //We will store the bit mask for the list of triggers.
-    //The list of triggers should be the same and in the same
-    //order, in order to read it later
-    Bool_t mIsGoodTrigger = false;
-    mCurrentTrigger = 0;
-    if(mTriggerIdCollection.empty()) {
+  //We will store the bit mask for the list of triggers.
+  //The list of triggers should be the same and in the same
+  //order, in order to read it later
+  Bool_t mIsGoodTrigger = false;
+  mCurrentTrigger = 0;
+  if(mTriggerIdCollection.empty()) {
+    mIsGoodTrigger = true;
+  }
+  else {
+    for (UInt_t iTrg = 0; iTrg < mTriggerIdCollection.size(); iTrg++) {
+      if(muEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) {
         mIsGoodTrigger = true;
-    }
-    else {
-        for (UInt_t iTrg = 0; iTrg < mTriggerIdCollection.size(); iTrg++) {
-            if(muEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) {
-                mIsGoodTrigger = true;
-                mCurrentTrigger |= (1<<iTrg);
-                break;
-            }
-        } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
-    }
+        mCurrentTrigger |= (1<<iTrg);
+        break;
+      }
+    } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
+  }
 
-    return mIsGoodTrigger;
+  return mIsGoodTrigger;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::AcceptPrimVtx(StThreeVectorF vtxPos, 
                                       Float_t vpdVz) {
-    Float_t mVtxX = vtxPos.x() - mPrimVtxXShift;
-    Float_t mVtxY = vtxPos.y() - mPrimVtxYShift;
-    Float_t vtxR = TMath::Sqrt(mVtxX*mVtxX + mVtxY*mVtxY);
-    Float_t vpdDiff = vtxPos.z() - vpdVz;
-
-    Bool_t mIsGoodVtx =  vtxPos.z() >= mPrimVtxZ[0] &&
-        vtxPos.z() <= mPrimVtxZ[1] &&
-        vtxR >= mPrimVtxR[0] &&
-        vtxR <= mPrimVtxR[1] &&
-        vpdDiff >= mPrimVtxVpdVzDiff[0] &&
-        vpdDiff <= mPrimVtxVpdVzDiff[1];
-
-    return mIsGoodVtx;
+  Float_t mVtxX = vtxPos.x() - mPrimVtxXShift;
+  Float_t mVtxY = vtxPos.y() - mPrimVtxYShift;
+  Float_t vtxR = TMath::Sqrt(mVtxX*mVtxX + mVtxY*mVtxY);
+  Float_t vpdDiff = vtxPos.z() - vpdVz;
+
+  Bool_t mIsGoodVtx =  vtxPos.z() >= mPrimVtxZ[0] &&
+    vtxPos.z() <= mPrimVtxZ[1] &&
+    vtxR >= mPrimVtxR[0] &&
+    vtxR <= mPrimVtxR[1] &&
+    vpdDiff >= mPrimVtxVpdVzDiff[0] &&
+    vpdDiff <= mPrimVtxVpdVzDiff[1];
+
+  return mIsGoodVtx;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) {
 
-    Bool_t mIsGoodTrack = trk->vertexIndex() == vtxInd &&
-        trk->momentum().mag() >= mTrackMomentum[0] && trk->momentum().mag() <= mTrackMomentum[1] &&
-        trk->type() == 1 &&
-        trk->eta() >= mTrackEta[0] && trk->eta() <= mTrackEta[1] &&
-        trk->nHits() >= mTrackNHits[0] && trk->nHits() <= mTrackNHits[1] &&
-        trk->nHitsFit() >= mTrackNHitsFit[0] && trk->nHitsFit() <= mTrackNHitsFit[1] &&
-        trk->flag() >= mTrackFlag[0] && trk->flag() <= mTrackFlag[1] &&
-        trk->dca(vtxInd).perp() >= mTrackDca[0] && trk->dca(vtxInd).perp() <= mTrackDca[1] &&
-        trk->dcaGlobal(vtxInd).perp() >= mTrackDcaGlobal[0] &&
-        trk->dcaGlobal(vtxInd).perp() <= mTrackDcaGlobal[1];
-
-    return mIsGoodTrack;
+  Bool_t mIsGoodTrack = trk->vertexIndex() == vtxInd &&
+    trk->momentum().mag() >= mTrackMomentum[0] && trk->momentum().mag() <= mTrackMomentum[1] &&
+    trk->type() == 1 &&
+    trk->eta() >= mTrackEta[0] && trk->eta() <= mTrackEta[1] &&
+    trk->nHits() >= mTrackNHits[0] && trk->nHits() <= mTrackNHits[1] &&
+    trk->nHitsFit() >= mTrackNHitsFit[0] && trk->nHitsFit() <= mTrackNHitsFit[1] &&
+    trk->flag() >= mTrackFlag[0] && trk->flag() <= mTrackFlag[1] &&
+    trk->dca(vtxInd).perp() >= mTrackDca[0] && trk->dca(vtxInd).perp() <= mTrackDca[1] &&
+    trk->dcaGlobal(vtxInd).perp() >= mTrackDcaGlobal[0] &&
+    trk->dcaGlobal(vtxInd).perp() <= mTrackDcaGlobal[1];
+
+  return mIsGoodTrack;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::AcceptRefMult2(StMuTrack *trk) {
 
-    Bool_t mIsGoodTrack = trk->flag()> 0 &&
-        trk->charge() != 0 &&
-        trk->nHitsFit() >= 10 &&
-        trk->dca().mag() < 3.0 &&
-        trk->momentum().mag() >= 1e-10 &&
-        TMath::Abs(trk->eta()) <= 1.0;
+  Bool_t mIsGoodTrack = trk->flag()> 0 &&
+    trk->charge() != 0 &&
+    trk->nHitsFit() >= 10 &&
+    trk->dca().mag() < 3.0 &&
+                       trk->momentum().mag() >= 1e-10 &&
+    TMath::Abs(trk->eta()) <= 1.0;
 
-    return mIsGoodTrack;
+  return mIsGoodTrack;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTofTrack(StMuTrack *trk) { //Only for globals
 
-    Bool_t mIsTofHit = trk->btofPidTraits().beta() > 0 &&
-        trk->btofPidTraits().timeOfFlight() > 0;
+  Bool_t mIsTofHit = trk->btofPidTraits().beta() > 0 &&
+    trk->btofPidTraits().timeOfFlight() > 0;
 
-    return mIsTofHit;
+  return mIsTofHit;
 }
 
 //_________________
 Int_t StFemtoDstMaker::Finish() {
 
-    mOutFile->cd();
-    mOutFile->Write();
-    //mTree->Write();
-    mOutFile->Close();
-
-    std::cout
-        << "*************************************" << "\n"
-        << "StFemtoDstMaker has been finished" << "\n"
-        << "\t nEventsPassed   : " << mNEventsPassed  << "\n"
-        << "\t nEventsProcessed: " << mNEventsIn << "\n"
-        << "Accept trigger fail = " << mAcceptTriggerFail << "\n"
-        << "RefMult fail = " << mRefMultFail << "\n"
-        << "Anti-pileup fail = " << mAntiPileupFail << "\n"
-        << "Ranking fail = " << mRankingFail << "\n"
-        << "VtxNull fail = " << mVtxNullFail << "\n"
-        << "NPrimVtx fail = " << mNPrimVtxFail << "\n"
-        << "AcceptPrimVtx fail = " << mAcceptPrimVtxFail << "\n"
-        << "*************************************" << "\n";
+  mOutFile->cd();
+  mOutFile->Write();
+  //mTree->Write();
+  mOutFile->Close();
+
+  std::cout
+    << "*************************************" << "\n"
+    << "StFemtoDstMaker has been finished" << "\n"
+    << "\t nEventsPassed   : " << mNEventsPassed  << "\n"
+    << "\t nEventsProcessed: " << mNEventsIn << "\n"
+    << "Accept trigger fail = " << mAcceptTriggerFail << "\n"
+    << "RefMult fail = " << mRefMultFail << "\n"
+    << "Anti-pileup fail = " << mAntiPileupFail << "\n"
+    << "Ranking fail = " << mRankingFail << "\n"
+    << "VtxNull fail = " << mVtxNullFail << "\n"
+    << "NPrimVtx fail = " << mNPrimVtxFail << "\n"
+    << "AcceptPrimVtx fail = " << mAcceptPrimVtxFail << "\n"
+    << "*************************************" << "\n";
 
-    return kStOk;
+  return kStOk;
 }
 
 //_________________
 void StFemtoDstMaker::CleanVariables() {
 
-    mEventIsGood = false;
+  mEventIsGood = false;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTpcPion(StMuTrack *pTrk) {
 
-    Bool_t isPion = false;
+  Bool_t isPion = false;
 
-    if( pTrk->momentum().mag() <= mPthresh && 
-       pTrk->nSigmaPion() >= mTpcPionNSigma[0] &&
-       pTrk->nSigmaPion() <= mTpcPionNSigma[1] ) {
-        isPion = true;
-    }
+  if( pTrk->momentum().mag() <= mPthresh && 
+      pTrk->nSigmaPion() >= mTpcPionNSigma[0] &&
+      pTrk->nSigmaPion() <= mTpcPionNSigma[1] ) {
+    isPion = true;
+  }
 
-    return isPion;
+  return isPion;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTpcKaon(StMuTrack *pTrk) {
 
-    Bool_t isKaon = false;
+  Bool_t isKaon = false;
 
-    if( pTrk->momentum().mag() <= mPthresh && 
-       pTrk->nSigmaKaon() >= mTpcKaonNSigma[0] &&
-       pTrk->nSigmaKaon() <= mTpcKaonNSigma[1] ) {
-        isKaon = true;
-    }
+  if( pTrk->momentum().mag() <= mPthresh && 
+      pTrk->nSigmaKaon() >= mTpcKaonNSigma[0] &&
+      pTrk->nSigmaKaon() <= mTpcKaonNSigma[1] ) {
+    isKaon = true;
+  }
 
-    return isKaon;
+  return isKaon;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTpcProton(StMuTrack *pTrk) {
 
-    Bool_t isProton = false;
+  Bool_t isProton = false;
 
-    if( pTrk->momentum().mag() <= mPthresh && 
-       pTrk->nSigmaProton() >= mTpcProtonNSigma[0] &&
-       pTrk->nSigmaProton() <= mTpcProtonNSigma[1] ) {
-        isProton = true;
-    }
+  if( pTrk->momentum().mag() <= mPthresh && 
+      pTrk->nSigmaProton() >= mTpcProtonNSigma[0] &&
+      pTrk->nSigmaProton() <= mTpcProtonNSigma[1] ) {
+    isProton = true;
+  }
 
-    return isProton;
+  return isProton;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTofPion(StMuTrack *gTrk) {
 
-    Bool_t isPion = false;
-    Float_t beta = -999.;
-    Float_t mSqr = -999.;
-    if( gTrk->btofPidTraits().beta() > 0 &&
-       gTrk->momentum().mag() > mPthresh ) {
-        beta = gTrk->btofPidTraits().beta();
-        mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
-
-        if(mSqr >= mTofPionMassSqr[0] &&
-           mSqr <= mTofPionMassSqr[1]) {
-            isPion = true;
-        }
+  Bool_t isPion = false;
+  Float_t beta = -999.;
+  Float_t mSqr = -999.;
+  if( gTrk->btofPidTraits().beta() > 0 &&
+      gTrk->momentum().mag() > mPthresh ) {
+    beta = gTrk->btofPidTraits().beta();
+    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
+
+    if(mSqr >= mTofPionMassSqr[0] &&
+       mSqr <= mTofPionMassSqr[1]) {
+      isPion = true;
     }
+  }
 
-    return isPion;
+  return isPion;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTofKaon(StMuTrack *gTrk) {
 
-    Bool_t isKaon = false;
-    Float_t beta = -999.;
-    Float_t mSqr = -999.;
-    if( gTrk->btofPidTraits().beta() > 0 &&
-       gTrk->momentum().mag() > mPthresh ) {
-        beta = gTrk->btofPidTraits().beta();
-        mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
-
-        if(mSqr >= mTofKaonMassSqr[0] &&
-           mSqr <= mTofKaonMassSqr[1]) {
-            isKaon = true;
-        }
+  Bool_t isKaon = false;
+  Float_t beta = -999.;
+  Float_t mSqr = -999.;
+  if( gTrk->btofPidTraits().beta() > 0 &&
+      gTrk->momentum().mag() > mPthresh ) {
+    beta = gTrk->btofPidTraits().beta();
+    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
+
+    if(mSqr >= mTofKaonMassSqr[0] &&
+       mSqr <= mTofKaonMassSqr[1]) {
+      isKaon = true;
     }
+  }
 
-    return isKaon;
+  return isKaon;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTofProton(StMuTrack *gTrk) {
 
-    Bool_t isProton = false;
-    Float_t beta = -999.;
-    Float_t mSqr = -999.;
-    if( gTrk->btofPidTraits().beta() > 0 &&
-       gTrk->momentum().mag() > 0.5 ) {
+  Bool_t isProton = false;
+  Float_t beta = -999.;
+  Float_t mSqr = -999.;
+  if( gTrk->btofPidTraits().beta() > 0 &&
+      gTrk->momentum().mag() > 0.5 ) {
 
-        //gTrk->momentum().mag() > mPthresh ) {
-        beta = gTrk->btofPidTraits().beta();
-        mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
+    //gTrk->momentum().mag() > mPthresh ) {
+    beta = gTrk->btofPidTraits().beta();
+    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
 
-        if(mSqr >= mTofProtonMassSqr[0] &&
-           mSqr <= mTofProtonMassSqr[1]) {
-            isProton = true;
-        }
+    if(mSqr >= mTofProtonMassSqr[0] &&
+       mSqr <= mTofProtonMassSqr[1]) {
+      isProton = true;
     }
+  }
 
-    return isProton;
-    }
+  return isProton;
+}
+
+//_________________
+void StFemtoDstMaker::SetTriggerId(const UInt_t &id) {
+  mTriggerIdCollection.push_back(id);
+}
 
-    //_________________
-    void StFemtoDstMaker::SetTriggerId(const UInt_t &id) {
-        mTriggerIdCollection.push_back(id);
+//_________________
+int StFemtoDstMaker::LoadEPCalibParam(int runid) {
+  if (!mCalibPath)   return 0;
+  
+  TFile *in = new TFile(Form("%s/OutMuAna_%d.root", mCalibPath, runid));
+  if (!in)   return 0;
+	
+  mQVMaker->LoadParameters(in);
+
+  in->Close();
+  delete in;
+  return 1;
+}
+
+//_________________
+int StFemtoDstMaker::GetCentralityBin(float refMult, int fVz) {
+  const int multMax = 1000;
+  const int mMultiplicityCut[2][17] = {{multMax, /* 0      */   // vpd-zdce-tac-protected
+                                        273,     /* 0   1  */   //-30<vz<-25
+                                        235,     /* 2   1  */
+                                        201,     /* 2   3  */
+                                        270,     /* 4   3  */
+                                        142,     /* 4   5  */
+                                        118,     /* 6   5  */
+                                        97,      /* 6   7  */
+                                        78,      /* 8   7  */
+                                        62,      /* 8   9  */
+                                        49,      /* 10  9  */
+                                        38,      /* 10  11 */
+                                        29,      /* 12  11 */
+                                        22,      /* 12  13 */
+                                        16,      /* 14  13 */
+                                        12,      /* 14  15 */
+                                        8},      /*     15 */
+                                       {multMax, // vz>-25
+                                        273, 
+                                        235, 
+                                        201, 
+                                        170, 
+                                        142, 
+                                        118, 
+                                        97, 
+                                        78, 
+                                        62, 
+                                        49, 
+                                        38, 
+                                        29, 
+                                        22, 
+                                        16, 
+                                        12, 
+                                        8}}; 
+
+  // NqvCent = 16 5% step up to 80%
+  int i_vz = fVz > 0 ? 1 : 0;
+  
+  for (int iCent = 1; iCent <= NqvCent; iCent++) {
+    if (mMultiplicityCut[i_vz][iCent] <= refMult &&
+        refMult < mMultiplicityCut[i_vz][iCent - 1]) {
+      
+      return iCent - 1;
     }
+  }
+  return -1;
+}
+
+//_________________
+void StFemtoDstMaker::CalcQvSMD(StMuEvent *mEv, QV fQv[NsubSMD], int fCent, int fVz) {
+
+  // Calculate raw flow vectors
+  float fqv[NsubSMD][4] = {{0}}; //east-west-combined, xy+weight
+  for (int iSub = 0; iSub < NsubSMD; iSub++) { // east-west
+
+    int nStrip;
+    for (int iXY = 0;iXY < 2; iXY++) {
+      if (iXY == 0)   nStrip = 8;  // vertical strips (x-direction)
+      else            nStrip = 9;  // horizontal strips (y-direction)
+
+      for (int iStrip = 1; iStrip < nStrip; iStrip++) {
+        // fqv[iSub][1-2] = position(cm) * adcThisTile
+        fqv[iSub][iXY]   += ZDCSMD_GetPosition(iSub, iXY, iStrip)*ZDCSMD(mEv, iSub, iXY, iStrip);
+        // fqv[iSub][3-4] = adcThisTile
+        fqv[iSub][iXY+2] += ZDCSMD(mEv, iSub, iXY, iStrip);
+      }
+
+      // normalization
+      if(fqv[iSub][iXY+2] > 0)
+        fqv[iSub][iXY] /= fqv[iSub][iXY+2];
+      else
+        fqv[iSub][iXY] = -9999;
+      
+      if (iSub == 0 && mEv->zdcTriggerDetector().adcSum(east) < 1)
+        fqv[iSub][iXY] = -9999;
+      if (iSub == 1 && mEv->zdcTriggerDetector().adcSum(west) < 1)
+        fqv[iSub][iXY] = -9999;
+      if (iSub == 2 && (fqv[0][iXY] < -9000 || fqv[1][iXY] < -9000))
+        fqv[iSub][iXY] = -9999;
+    } // iXY
+  } // iSub
+  
+  // combined
+  for (int iXY = 0; iXY < 2; iXY++) {
+    fqv[2][iXY] = fqv[0][iXY] - fqv[1][iXY];
+    if (fqv[0][iXY] == -9999 || fqv[1][iXY] == -9999)
+      fqv[2][iXY] = -9999;
+  }
+
+  
+  // Set raw Q-vectors and Calibration 
+  // -----------------------------------------------
+  int fth   = 0; // fth = fOrd in QV class
+  float fQw = 1.0;
+  for (int iSub = 0; iSub < NsubSMD; iSub++) {
+    if (fqv[iSub][0] == -9999 || fqv[iSub][1] == -9999 )   fQw = -1;
+    fQv[iSub].set_Qv(fqv[iSub][0], fqv[iSub][1], fQw, dZDC, iSub, fth, fCent, fVz);
+
+    mQVMaker->doCalibration(&(fQv[iSub]));
+  } // iSub
+}
 
+//_________________
+Float_t StFemtoDstMaker::ZDCSMD(StMuEvent *mEv, int eastwest, int verthori, int strip) {
+  
+  // return adc for each tile
+  return mEv->zdcTriggerDetector().zdcSmd((StBeamDirection)eastwest, verthori, strip);
+}
+
+//_________________
+Float_t StFemtoDstMaker::ZDCSMD_GetPosition(int eastwest, int verthori, int strip) {
+
+  // Get position of each slat;strip starts from 1
+  Float_t zdcsmd_x[7] = {0.5, 2, 3.5, 5, 6.5, 8, 9.5};
+  Float_t zdcsmd_y[8] = {1.25, 3.25, 5.25, 7.25, 9.25, 11.25, 13.25, 15.25};
+
+  if (eastwest == 0 && verthori == 0) { return zdcsmd_x[strip-1] - mZDCSMDCenterex; }
+  // east && vertical strips (x-direction)
+  if (eastwest == 1 && verthori == 0) { return mZDCSMDCenterwx - zdcsmd_x[strip-1]; }
+  // west && vertical strips (x-direction)
+  if (eastwest == 0 && verthori == 1) { return zdcsmd_y[strip-1]/sqrt(2.) - mZDCSMDCenterey; }
+  // east && horisontal strips (y-direction)
+  if (eastwest == 1 && verthori == 1) { return zdcsmd_y[strip-1]/sqrt(2.) - mZDCSMDCenterwy; }
+  // west && horisontal strips (y-direction)
+  
+  return 0;
+}

+ 22 - 1
StRoot/StFemtoDstMaker/StFemtoDstMaker.h

@@ -12,6 +12,9 @@
 #include "StPhysicalHelixD.hh"
 #include "StRefMultCorr/StRefMultCorr.h"
 
+#include "QVMaker/QV.h"
+#include "QVMaker/QVMaker.h"
+
 #include "StFemtoEvent.h"
 #include "StFemtoTrack.h"
 
@@ -55,6 +58,7 @@ class StMuTrack;
 class StFemtoDstMaker : public StMaker {
  public:
   StFemtoDstMaker(StMuDstMaker *muMaker,
+                  const char *calibPath,
                   const Char_t *oFileName);
   ~StFemtoDstMaker();
 
@@ -116,6 +120,23 @@ class StFemtoDstMaker : public StMaker {
 
   Int_t  mNBytes;                      // number of writed bytes to the mTree
   vector<UInt_t> mTriggerIdCollection;
+
+  //
+  // Reaction plane
+  //
+  QVMaker *mQVMaker;
+  double   mZDCSMDCenterex = 4.72466; //
+  double   mZDCSMDCenterey = 5.53629; // This variables
+  double   mZDCSMDCenterwx = 4.39604; //  is from Takafumi codes
+  double   mZDCSMDCenterwy = 5.19968; //
+  
+  const char *mCalibPath; // path to the calibration ROOT files
+
+  int LoadEPCalibParam(int runid); // load calibration
+  int GetCentralityBin(float refMult, int fVz);
+  void CalcQvSMD(StMuEvent *mEv, QV fQv[NsubSMD], int fCent, int fVz);
+  Float_t ZDCSMD(StMuEvent *mEv, int eastwest, int verthori, int strip);
+  Float_t ZDCSMD_GetPosition(int eastwest, int verthori, int strip);
   //
   // Transverse sphericity
   //
@@ -211,7 +232,7 @@ class StFemtoDstMaker : public StMaker {
   Bool_t IsTpcProton(StMuTrack *pTrk);
   Bool_t IsTofProton(StMuTrack *gTrk);
 
-  ClassDef(StFemtoDstMaker, 3)
+  ClassDef(StFemtoDstMaker, 4)
 };
 
 #endif

+ 9 - 1
StRoot/StFemtoDstMaker/StFemtoDstQAMaker.cxx

@@ -133,6 +133,12 @@ Int_t StFemtoDstQAMaker::Init() {
                                        mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
     hTriggerId = new TH1I("hTriggerId", ";trigger bit;#", 32, 0, 32);
     hSphericity = new TH1F("hSphericity", ";sphericity;#", 120, -0.1, 1.1);
+    hRPeast = new TH1F("hRPeast", ";#Psi_{east};#frac{dN}{d#Psi_{east}}",
+                       200, -4, 4.);
+    hRPwest = new TH1F("hRPwest", ";#Psi_{west};#frac{dN}{d#Psi_{west}}",
+                       200, -4, 4.);
+    hRPdiff = new TH1F("hRPdiff", ";#Psi_{east} - #Psi_{west};#frac{dN}{d(#Psi_{east} - #Psi_{west}}",
+                       700, -7, 7.);    
     // Track histograms
     hNSigmaPionVsMom = new TH2F("hNSigmaPionVsMom", "N#sigma(#pi) vs momenta;p (GeV/c);N#sigma(pi)",
                                 mMomBins,mMomLo,mMomHi,
@@ -399,7 +405,9 @@ Int_t StFemtoDstQAMaker::Make() {
         hRefMult2VsRunNumber->Fill(femtoEvent->GetRunId(),
                                    femtoEvent->GetRefMult2());
         hSphericity->Fill(femtoEvent->GetSphericity());
-
+        hRPeast->Fill(femtoEvent->GetRPeast());
+        hRPwest->Fill(femtoEvent->GetRPwest());
+        hRPdiff->Fill(femtoEvent->GetRPeast() - femtoEvent->GetRPwest());
         TClonesArray *lTracks = femtoEvent->GetFemtoTracks();
 
         if (lTracks) {

+ 3 - 0
StRoot/StFemtoDstMaker/StFemtoDstQAMaker.h

@@ -64,6 +64,9 @@ class StFemtoDstQAMaker : public StMaker {
   TH2F *hPrimVertXvsY;
   TH1F *hPrimVertVpdVzDiff;
   TH1I *hTriggerId;
+  TH1F *hRPeast;
+  TH1F *hRPwest;
+  TH1F *hRPdiff;
   //Here is the way of how to look on the quantities over the run
   TProfile *hRefMultVsRunNumber;
   TProfile *hPrimTrackNumVsRunNumber;

+ 31 - 26
StRoot/StFemtoDstMaker/StFemtoEvent.h

@@ -46,46 +46,51 @@ class StFemtoEvent : public TObject {
   void SetTriggerId(UInt_t id)        { mTriggerId = id; }
   void SetPrimaryVertexRanking(Float_t);
   void SetSphericity(Float_t sph)     { mSphericity = (char)round(sph*100.); }
+  void SetRPeast(Float_t rp)              { mRPeast = rp; }
+  void SetRPwest(Float_t rp)              { mRPwest = rp; }  
   StFemtoTrack* AddFemtoTrack();
   void Clear(Option_t* option = "C");
 
   //Getters
-  Int_t    GetEventId() const          { return mEventId; }
-  Int_t    GetRunId() const            { return mRunId; }
-  Bool_t   GetCollisionType() const    { return mCollisionType; }
-  UInt_t   GetRefMult() const          { return (UInt_t)mRefMult; }
-  UInt_t   GetRefMult2() const         { return (UInt_t)mRefMult2;}
-  Float_t  GetRefMultCorr() const      { return (Float_t)mRefMultCorr * 0.1; }
-  Float_t  GetRefMultCorrWeight() const{ return (Float_t)mRefMultCorrWeight * 0.0001; }
+  Int_t    GetEventId() const               { return mEventId; }
+  Int_t    GetRunId() const                 { return mRunId; }
+  Bool_t   GetCollisionType() const         { return mCollisionType; }
+  UInt_t   GetRefMult() const               { return (UInt_t)mRefMult; }
+  UInt_t   GetRefMult2() const              { return (UInt_t)mRefMult2;}
+  Float_t  GetRefMultCorr() const           { return (Float_t)mRefMultCorr * 0.1; }
+  Float_t  GetRefMultCorrWeight() const     { return (Float_t)mRefMultCorrWeight * 0.0001; }
   UShort_t GetCent9() const;
-  UShort_t GetCent16() const           { return (UShort_t)mCent16; }
-  UShort_t GetRefMultPos() const       { return mRefMultPos; }
-  UShort_t GetRefMultNeg() const       { return (mRefMult - mRefMultPos); }
-  UShort_t GetRefMult2Pos() const       { return mRefMult2Pos; }
-  UShort_t GetRefMult2Neg() const       { return (mRefMult2 - mRefMult2Pos); }
-  Float_t  GetZDCe() const             { return (Float_t)30000; }
-  Float_t  GetZDCw() const             { return (Float_t)30000; }
-  UShort_t GetNumberOfBTofHit() const  { return mNumberOfBTofHit; }
+  UShort_t GetCent16() const                { return (UShort_t)mCent16; }
+  UShort_t GetRefMultPos() const            { return mRefMultPos; }
+  UShort_t GetRefMultNeg() const            { return (mRefMult - mRefMultPos); }
+  UShort_t GetRefMult2Pos() const           { return mRefMult2Pos; }
+  UShort_t GetRefMult2Neg() const           { return (mRefMult2 - mRefMult2Pos); }
+  Float_t  GetZDCe() const                  { return (Float_t)30000; }
+  Float_t  GetZDCw() const                  { return (Float_t)30000; }
+  UShort_t GetNumberOfBTofHit() const       { return mNumberOfBTofHit; }
   UShort_t GetNumberOfPrimaryTracks() const { return mNumberOfPrimaryTracks; }
   UShort_t GetNumberOfGlobalTracks() const  { return mNumberOfGlobalTracks; }
-  Float_t  GetMagneticField() const    { return (Float_t)mMagField*0.001; }
+  Float_t  GetMagneticField() const         { return (Float_t)mMagField*0.001; }
   TVector3 GetPrimaryVertexPosition() const;
-  Float_t  GetVertexPositionX() const  { return mVertexPositionX; }
-  Float_t  GetVertexPositionY() const  { return mVertexPositionY; }
-  Float_t  GetVertexPositionZ() const  { return mVertexPositionZ; }
+  Float_t  GetVertexPositionX() const       { return mVertexPositionX; }
+  Float_t  GetVertexPositionY() const       { return mVertexPositionY; }
+  Float_t  GetVertexPositionZ() const       { return mVertexPositionZ; }
   Float_t  GetVertexPositionR() const;
-  Float_t  GetVpdVz() const            { return mVpdVz; }
-  UInt_t   GetTriggerId() const        { return mTriggerId; }
-  Float_t  GetPrimaryVertexRanking() const {return (Float_t)mPrimaryVertexRanking; }
-  TClonesArray *GetFemtoTracks() const { return mFemtoTracks; }
-  Int_t    GetNFemtoTracks() const { return (Int_t)mNFemtoTracks; }
-  Float_t  GetSphericity() const { return ((float)mSphericity)/100.; }
+  Float_t  GetVpdVz() const                 { return mVpdVz; }
+  UInt_t   GetTriggerId() const             { return mTriggerId; }
+  Float_t  GetPrimaryVertexRanking() const  { return (Float_t)mPrimaryVertexRanking; }
+  TClonesArray *GetFemtoTracks() const      { return mFemtoTracks; }
+  Int_t    GetNFemtoTracks() const          { return (Int_t)mNFemtoTracks; }
+  Float_t  GetSphericity() const            { return ((float)mSphericity)/100.; }
+  Float_t  GetRPeast() const                    { return mRPeast; }
+  Float_t  GetRPwest() const                    { return mRPwest; }
 
  private:
 
   Int_t    mEventId;
   Int_t    mRunId;
   Char_t   mSphericity;          //transverse sphericity
+  Float_t  mRPeast, mRPwest;     // reaction plane
   Bool_t   mCollisionType;       //0-pp, 1-AA(from those where RefMultCorr is defined)
   UShort_t mRefMult;             //mRefMult
   UShort_t mRefMultCorr;         //mRefMultCorr * 10
@@ -108,7 +113,7 @@ class StFemtoEvent : public TObject {
   UShort_t mNFemtoTracks;
   TClonesArray *mFemtoTracks;
 
-  ClassDef(StFemtoEvent, 101)
+  ClassDef(StFemtoEvent, 102)
 };
 
 #endif

+ 2 - 0
StRoot/StFemtoDstMaker/StHbtFemtoDstReader.cxx

@@ -231,6 +231,8 @@ StHbtEvent *StHbtFemtoDstReader::Read() {
         mHbtEvent->SetRefMultCorr( mFemtoEvent->GetRefMultCorr() );
         mHbtEvent->SetSphericity( mFemtoEvent->GetSphericity() );
         mHbtEvent->SetTofMult( mFemtoEvent->GetNumberOfBTofHit() );
+        mHbtEvent->SetPsiRP1(0, mFemtoEvent->GetRPeast() );
+        mHbtEvent->SetPsiRP1(1, mFemtoEvent->GetRPwest() );
 
         TClonesArray *mFemtoTrackArr = mFemtoEvent->GetFemtoTracks();
 

+ 40 - 20
StRoot/StFemtoDstMaker/StHbtO97DstReader.cxx

@@ -15,6 +15,7 @@ StHbtO97DstReader::StHbtO97DstReader(const char *dir, const char *fileName,
 
     mTotalTracks = 0.;
 
+    mIsaJet     = 0; // default is PDG codes
     mTree       = 0;
     mTChain     = 0;
     mO97Event   = 0;
@@ -32,6 +33,15 @@ int StHbtO97DstReader::InitRead(string dir, string fileName,
 
     std::cout << "[StHbtO97DstReader]: ---=== StHbtO97DstReader initialization ===---" << std::endl;
 
+    //
+    // Load PDG database
+    //
+    pdgDb = new TDatabasePDG;
+    if (pdgDb) 
+      pdgDb->ReadPDGTable("StRoot/StFemtoDstMaker/pdg_table.txt");
+    else
+      return 0;
+
     int numFiles = 0;
     if (!fileName.empty()) {
 
@@ -186,8 +196,8 @@ StHbtEvent *StHbtO97DstReader::Read() {
         mHbtEvent = new StHbtEvent();
 
         int nTracks = mO97Event->GetNTracks();
-        int nPosTracks = 0;
-        int nNegTracks = 0;
+        int refMultPos = 0;
+        int refMultNeg = 0;
         StThreeVectorF vPos( 0., 0., 0. );
 
         TClonesArray *mO97TrackArr = mO97Event->GetTracks();
@@ -196,7 +206,7 @@ StHbtEvent *StHbtO97DstReader::Read() {
             mTotalTracks += (float)nTracks;
             if (mDebug) {
                 std::cout << "[StHbtO97DstReader]: nTracks = "
-                    << nTracks << std::endl;
+                          << nTracks << std::endl;
             }
 
             for (int i = 0; i < nTracks; i++) {
@@ -206,25 +216,36 @@ StHbtEvent *StHbtO97DstReader::Read() {
                 StHbtTrack *mHbtTrack = new StHbtTrack();
 
                 // Set charge, mass and momentum
-                char charge = mO97Track->GetPdgId()/abs(mO97Track->GetPdgId());
-                if (charge > 0) {
-                    nPosTracks++;
-                }
-                else {
-                    nNegTracks++;
-                }
+                char charge = mO97Track->GetCharge(pdgDb, mIsaJet);
 
                 float trkMass = mO97Track->GetMass();
                 StThreeVectorF trkMom(mO97Track->GetPx(), 
                                       mO97Track->GetPy(), 
                                       mO97Track->GetPz());
+                // Calculate DCA in cm
+                StThreeVectorF origin(mO97Track->GetXfr()*1e-13,
+                                      mO97Track->GetYfr()*1e-13,
+                                      mO97Track->GetZfr()*1e-13);
+                float dcaXY = origin.perp();
+                float dcaZ = origin.z();
 
+                // Calculate reference multiplicity
+                if (charge != 0 &&
+                    fabs(trkMom.pseudoRapidity()) < 0.5 &&
+                    trkMom.mag() >= 1e-10 &&
+                    origin.mag() < 3.)
+                {
+                  if (charge > 0)   refMultPos++;
+                  else              refMultNeg++;
+                }
+                
                 // Set nSigma and PID probability
                 float trkNSigmaElectron,  trkNSigmaPion;
                 float trkNSigmaKaon,      trkNSigmaProton;
                 float trkPidProbElectron, trkPidProbPion;
-                float trkPidProbKaon,     trkPidProbProton;	
-                switch (abs(mO97Track->GetPdgId())) {
+                float trkPidProbKaon,     trkPidProbProton;
+                int pdg = mIsaJet ? mO97Track->GetPdgId(pdgDb) : mO97Track->GetPdgId();
+                switch (abs(pdg)) {
                     case 211: // pion
                         trkNSigmaElectron = -99.;
                         trkNSigmaPion     = 0.;
@@ -290,11 +311,6 @@ StHbtEvent *StHbtO97DstReader::Read() {
                 mHbtTrack->SetdEdx( dedxMean(trkMass, trkMom.mag()) );
 
                 // Set DCA in cm
-                StThreeVectorF origin( mO97Track->GetXfr()*1e-13,
-                                      mO97Track->GetYfr()*1e-13,
-                                      mO97Track->GetZfr()*1e-13 );
-                float dcaXY = origin.perp();
-                float dcaZ = origin.z();
                 mHbtTrack->SetDCAxyGlobal( dcaXY );
                 mHbtTrack->SetDCAzGlobal( dcaZ );
                 mHbtTrack->SetChi2(1.);
@@ -322,6 +338,9 @@ StHbtEvent *StHbtO97DstReader::Read() {
                 mHbtTrack->SetFreezeOutZ( mO97Track->GetZfr() );
                 mHbtTrack->SetFreezeOutT( mO97Track->GetTfr() );
 
+                mHbtTrack->SetPdgCode( mO97Track->GetPdgId() );
+                mHbtTrack->SetSimTrack();
+
                 mHbtEvent->TrackCollection()->push_back( mHbtTrack );
             } // for (int i = 0; i < nTracks; i++)
         } // if (mO97TrackArr)
@@ -336,9 +355,10 @@ StHbtEvent *StHbtO97DstReader::Read() {
         mHbtEvent->SetReactionPlane( 0., 0 );
         mHbtEvent->SetReactionPlaneError( 0., 0 );
         mHbtEvent->SetReactionPlaneSubEventDifference( 0., 0 );
-        mHbtEvent->SetUncorrectedNumberOfPositivePrimaries( nPosTracks );
-        mHbtEvent->SetUncorrectedNumberOfNegativePrimaries( nNegTracks );
-        mHbtEvent->SetUncorrectedNumberOfPrimaries( nTracks );
+        mHbtEvent->SetUncorrectedNumberOfPositivePrimaries( refMultPos );
+        mHbtEvent->SetUncorrectedNumberOfNegativePrimaries( refMultNeg );
+        mHbtEvent->SetUncorrectedNumberOfPrimaries( refMultPos + refMultNeg );
+        mHbtEvent->SetTofMult( refMultPos + refMultNeg );
         mHbtEvent->SetPrimVertPos( vPos );
         mHbtEvent->SetMagneticField( PYTHIA_RBFIELD );
         mHbtEvent->SetVpdVz( 0. );

+ 11 - 7
StRoot/StFemtoDstMaker/StHbtO97DstReader.h

@@ -18,6 +18,7 @@
 #include "StChain.h"
 #include <string>
 #include <stdlib.h>
+#include "TDatabasePDG.h"
 
 class StFlowEvent;
 class StIOMaker;
@@ -26,12 +27,14 @@ class TFile;
 
 //_________________
 class StHbtO97DstReader : public StHbtEventReader {
- private:
+private:
   string mDir;
   string mFileName;
   string mFilter;
 
+  TDatabasePDG *pdgDb;  
   bool          mDebug;
+  bool          mIsaJet;      // default is the PDG codes
   unsigned int  mEventIndex;
   unsigned int  mNEvents;
   TChain       *mTChain;
@@ -42,23 +45,24 @@ class StHbtO97DstReader : public StHbtEventReader {
   int mMaxFiles;
 
   int InitRead(string dir, string fileName,
-	       string filter, int maxFiles);
+               string filter, int maxFiles);
   int FillChain(TChain *chain, const char *fileName, 
-		int maxFiles);
+                int maxFiles);
   int FillChain(TChain *chain, char *dir,
-		const char *filter, int maxFiles);
+                const char *filter, int maxFiles);
   void UninitRead();
   float dedxMean(float mass, float mom);
   StHbtEvent *Read();
- public:
+public:
   StHbtO97DstReader(const char *dir, const char *fileName,
-		    const char *filter = ".", int maxFiles = 1e9, bool debug = false);
+                    const char *filter = ".", int maxFiles = 1e9, bool debug = false);
   StHbtEvent *ReturnHbtEvent();
   void PrintTotalTracks();
   StHbtString Report();
   int GetNEvents();
+  void SetIsaJetMode(bool mode) { mIsaJet = mode; }
 
-  ClassDef(StHbtO97DstReader, 1)
+  ClassDef(StHbtO97DstReader, 2)
 };
 
 #endif

+ 34 - 21
StRoot/StFemtoDstMaker/StO97DstMaker.cxx

@@ -17,8 +17,14 @@ StO97DstMaker::StO97DstMaker(const Char_t *iFileName,
     TDatime d;
 
     unsigned int num;
-    sscanf(oFileName, "%*[^_]_%d.f19", &num);
-    mOutFileName = Form("%s_%d%d%d", oFileName, d.GetYear(), d.GetMonth(), num);
+    TString ifname(iFileName);
+    int dotPos = ifname.Last('.');
+    ifname[dotPos] = '\0';
+    num = atoi(ifname.Data() + ifname.Last('_') + 1);
+    ifname[dotPos] = '.';
+    mOutFileName = Form("%s_%04d%02d%02d_%d.root", oFileName,
+                        d.GetYear(), d.GetMonth(), d.GetDay(), num);
+    std::cout << "[StO97DstMaker]: output file name " << mOutFileName << std::endl;
 
     mInFileName  = iFileName;
     mTree = 0;
@@ -28,7 +34,8 @@ StO97DstMaker::StO97DstMaker(const Char_t *iFileName,
     mNEvents = 0;
     mStop = false;
 
-    matrix = new TMatrixTSym<double>(2); // matrix for calculation of transverse sphericity
+    matrix  = new TMatrixTSym<double>(2); // matrix for calculation of transverse sphericity
+    matrix2 = new TMatrixTSym<double>(2); // matrix for calculation of transverse sphericity 2
 
     std::cout << "[OK]" << std::endl;
 }
@@ -48,11 +55,13 @@ Int_t StO97DstMaker::Init() {
     // Input file
     //
     std::cout << "[StO97DstMaker] Openning an input file... ";
-    if (mInFile)
+    if (mInFile) {
         fclose(mInFile);
+    }
     mInFile = fopen(mInFileName, "r");
-    if (mInFile)
+    if (mInFile) {
         std::cout << "[OK]" << std::endl;
+    }
     else {
         std::cout << "[FAIL]" << std::endl;
         return kStFatal;
@@ -65,10 +74,12 @@ Int_t StO97DstMaker::Init() {
         mOutFile->Close();
         mOutFile->Open(mOutFileName, "RECREATE");
     }
-    else
+    else {
         mOutFile = new TFile(mOutFileName, "RECREATE");
-    if (mOutFile)
+    }
+    if (mOutFile) {
         std::cout << "[OK]" << std::endl;
+    }
     else {
         std::cout << "[FAIL]" << std::endl;
         return kStFatal;
@@ -78,8 +89,9 @@ Int_t StO97DstMaker::Init() {
     // TTree setup
     //
     std::cout << "[StO97DstMaker] Creating TTree... ";
-    if (mTree)
+    if (mTree) {
         delete mTree;
+    }
     mTree = new TTree("StO97Dst", "StO97Dst");
     if (!mTree) {
         std::cout << "[FAIL]" << std::endl;
@@ -112,9 +124,11 @@ Int_t StO97DstMaker::Make() {
     // Read event
     //
     mNEvents++;
+    mStop = !mStop;
     mStop = ReadEvent();
-    if (mStop)
+    if (mStop) {
         return kStOk;
+    }
     //
     // Read tracks
     //
@@ -143,8 +157,9 @@ Int_t StO97DstMaker::Finish() {
         mInFile = 0;
     }
     if (mTree) {
-        if (mOutFile)
+        if (mOutFile) {
             mOutFile->Write();
+        }
         delete mTree;
         mTree = 0;
     }
@@ -161,31 +176,29 @@ Int_t StO97DstMaker::Finish() {
 //_________________
 Bool_t StO97DstMaker::ReadEvent() {
     Int_t eventNumber;
-    Int_t nTracks;
     Float_t impPar;
     Float_t rot;
     Char_t buf[128];
 
-
     while (1) {
         Int_t ret = fscanf(mInFile, " %i %i %f %f",
-                           &eventNumber, &nTracks, &impPar, &rot);
-        if (ret <= 0)               // if can't read a line
+                           &eventNumber, &mNTracks, &impPar, &rot);
+        if (ret <= 0) {               // if can't read a line
             fgets(buf, 127, mInFile); // skip it
+        }
         else if (ret == 4) { // 4 - total number of event parameters
             mEvent->Clear();
             mEvent->SetEventNumber(eventNumber);
-            mEvent->SetNTracks(nTracks);
             mEvent->SetImpactPar(impPar);
             mEvent->SetEventPlaneRot(rot);
             return false;
         }
-        if (feof(mInFile)) // if EOF then return stop flag
+        if (feof(mInFile)) { // if EOF then return stop flag
             return true;
+        }
     }
 }
 
-
 //_________________
 Bool_t StO97DstMaker::ReadTracks() {
     Int_t id;
@@ -193,7 +206,6 @@ Bool_t StO97DstMaker::ReadTracks() {
     Float_t px, py, pz;
     Float_t energy, mass;
     Float_t xFr, yFr, zFr, tFr;
-    Int_t nTracks = mEvent->GetNTracks();
     TLorentzVector particle;
 
     Float_t pt_full  = 0.0;
@@ -204,12 +216,13 @@ Bool_t StO97DstMaker::ReadTracks() {
     matrix->Zero();
     matrix2->Zero();
 
-    for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+    for (Int_t iTrack = 0; iTrack < mNTracks; iTrack++) {
         Int_t ret = fscanf(mInFile, " %i %i %f %f %f %f %f %f %f %f %f",
                            &id, &pdgId, &px, &py, &pz,
                            &energy, &mass, &xFr, &yFr, &zFr, &tFr);
-        if (ret != 11)
+        if (ret != 11) {
             return true;
+        }
         else {
             StO97Track *track = mEvent->AddOscar97Track();
             track->SetId(id);
@@ -230,10 +243,10 @@ Bool_t StO97DstMaker::ReadTracks() {
              */
             // setup particle vector
             particle.SetPxPyPzE(px, py, pz, energy);
-            Float_t eta = particle.PseudoRapidity();
             Float_t pt  = particle.Pt();
 
             if (pt > 0.15) {
+                Float_t eta = particle.PseudoRapidity(); // it heals pt = 0 case
                 if (TMath::Abs(eta) < 0.5)
                 {
                     (*matrix)(0, 0) += px*px/pt;

+ 1 - 0
StRoot/StFemtoDstMaker/StO97DstMaker.h

@@ -54,6 +54,7 @@ class StO97DstMaker : public StMaker {
   Bool_t mStop;
   Int_t  mNBytes;
   Int_t  mNEvents;
+  Int_t  mNTracks;
   char   mBuf[128]; // buffer for general read purposes
 
   StO97Event *mEvent;

+ 305 - 252
StRoot/StFemtoDstMaker/StO97DstQAMaker.cxx

@@ -4,303 +4,356 @@ ClassImp(StO97DstQAMaker)
 
 //_________________
 StO97DstQAMaker::StO97DstQAMaker(const char *dirName,
-				 const char *fileName,
-				 const char *filter,
-				 int        maxFiles) {
-   mDir = string(dirName);
-   mFileName = string(fileName);
-   mFilter = string(filter);
-
-   mTree     = 0;
-   mChain    = 0;
-   mO97Event = 0;
-   mEventId  = 0;
+                                 const char *fileName,
+                                 const char *filter,
+                                 int        maxFiles) {
+    mDir = string(dirName);
+    mFileName = string(fileName);
+    mFilter = string(filter);
+
+    mTree     = 0;
+    mChain    = 0;
+    mO97Event = 0;
+    mEventId  = 0;
 }
- 
+
 
 //_________________
 Int_t StO97DstQAMaker::Init() {
-  std::cout << "[StO97DstQAMaker]: Initialization... ";
-  InitRead(mDir, mFileName, mFilter, mMaxFiles);
-  mNEvents = (unsigned int)mChain->GetEntries();
-  mOutFile = new TFile(mOutFileName, "RECREATE");
-  if (!mOutFile) {
-    std::cout << "FAIL" << std::endl;
-    return kStErr;
-  }
-  else
-    std::cout << "OK" << std::endl;
-  //
-  // Load PDG database
-  //
-  pdgDb = new TDatabasePDG;
-  pdgDb->ReadPDGTable("/star/u/pusheax/pdg_table.txt");
-  //
-  // Creating event histograms
-  //
-  hRefMultP90  = new TH1F("hRefMultP90", "reference multiplicity (greg) P_{t} > 0.09 GeV/c;RefMult;#", 20, 0., 20.);
-  hRefMultP100 = new TH1F("hRefMultP100", "reference multiplicity (greg) P_{t} > 0.1 GeV/c;RefMult;#", 20, 0., 20.);
-  hRefMultP120 = new TH1F("hRefMultP120", "reference multiplicity (greg) P_{t} > 0.12 GeV/c;RefMult;#", 20, 0., 20.);
-  hRefMultP150 = new TH1F("hRefMultP150", "reference multiplicity (greg) P_{t} > 0.15 GeV/c;RefMult;#", 20, 0., 20.);
-  hImpPar      = new TH1F("hImpPar", "impact parameter;imp (fm/c);#", 100, 0., 20.);
-  hImpParVsRefMNik = new TH2F("hImpParVsRefMNik",
-			      "impact parameter vs reference multiplicity (nik)",
-			      100, 0., 13.,      // impact parameter
-			      2000, 0., 20000.); // reference multiplicity
-  hImpParVsRefMGreg = new TH2F("hImpParVsRefMGreg",
-			       "impact parameter vs reference multiplicity (greg)",
-			       100, 0., 13.,      // impact parameter
-			       2000, 0., 20000.); // reference multiplicity
-  //
-  // Creating track histograms
-  //
-  hVx = new TH1F("hVx", "x freeze-out position;V_{x} (fm);#", 100, -2., 2.);
-  hVy = new TH1F("hVy", "y freeze-out position;V_{y} (fm);#", 100, -2., 2.);
-  hVz = new TH1F("hVz", "z freeze-out position;V_{z} (fm);#", 100, -2., 2.);
-  hPx = new TH1F("hPx", "single particle p_{x};p_{x} (GeV/c);#", 200, -2., 2.);
-  hPy = new TH1F("hPy", "single particle p_{y};p_{y} (GeV/c);#", 200, -2., 2.);
-  hPz = new TH1F("hPz", "single particle p_{z};p_{z} (GeV/c);#", 200, -2., 2.);
-  hP  = new TH1F("hP", "single particle p;p (GeV/c);#", 400, 0., 2.);
-  hPt = new TH1F("hPt", "single particle p_{t};p_{t} (GeV/c);#", 400, 0., 2.);  
-  hEta = new TH1F("hEta", ";#eta;#", 200, -10., 10.);
-  hEtaInt = new TH1F("hEtaInt", ";#eta;#", 200, -10., 10.);  
-  hMSqrVsPt = new TH2F("hMSqrVsPt",
-		       "m^{2} vs p__{t}:p_{t}*charge (GeV/c):m^{2} (GeV^{2})",
-		       200, -2., 2., // Pt
-		       100, 0., 2.); // msqr
-
-  return kStOk;
+    std::cout << "[StO97DstQAMaker]: Initialization... ";
+    InitRead(mDir, mFileName, mFilter, mMaxFiles);
+    mNEvents = (unsigned int)mChain->GetEntries();
+    mOutFile = new TFile(mOutFileName, "RECREATE");
+    if (!mOutFile) {
+        std::cout << "FAIL" << std::endl;
+        return kStErr;
+    }
+    else
+        std::cout << "OK" << std::endl;
+    //
+    // Load PDG database
+    //
+    pdgDb = new TDatabasePDG;
+    pdgDb->ReadPDGTable("StRoot/StFemtoDstMaker/pdg_table.txt");
+    //
+    // Creating event histograms
+    //
+		int refMultBins = 500;
+		float refMultLo = 0.;
+		float refMultHi = 500;
+    hRefMultP90  = new TH1F("hRefMultP90", "reference multiplicity (greg) p_{#perp} > 0.09 GeV/c;RefMult;#",  refMultBins, refMultLo, refMultHi);
+    hRefMultP100 = new TH1F("hRefMultP100", "reference multiplicity (greg) p_{#perp} > 0.1 GeV/c;RefMult;#",  refMultBins, refMultLo, refMultHi);
+    hRefMultP120 = new TH1F("hRefMultP120", "reference multiplicity (greg) p_{#perp} > 0.12 GeV/c;RefMult;#", refMultBins, refMultLo, refMultHi);
+    hRefMultP150 = new TH1F("hRefMultP150", "reference multiplicity (greg) p_{#perp} > 0.15 GeV/c;RefMult;#", refMultBins, refMultLo, refMultHi);
+    hImpPar      = new TH1F("hImpPar", "impact parameter;imp (fm/c);#", 100, 0., 20.);
+    hSph         = new TH1F("hSph",
+                            "Transverse sphericity (|#eta| < 0.5, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
+                            100, 0., 1.);
+    hSph         = new TH1F("hSph",
+                            "Transverse sphericity (|#eta| < 0.5, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
+                            100, 0., 1.);
+    hSph2        = new TH1F("hSph2",
+                            "Transverse sphericity (|#eta| < 1.0, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
+                            100, 0., 1.);
+		hImpVsRef    = new TH2F("hImpVsRef", "p_{#perp} > 0.15 [GeV/c] && |#eta| < 0.5;b [fm];RefMult", 100, 0., 20., refMultBins, refMultLo, refMultHi);
+
+    //
+    // Creating track histograms
+    //
+    hVx = new TH1F("hVx", "x freeze-out position;V_{x} (fm);#", 500, -50., 50.);
+    hVy = new TH1F("hVy", "y freeze-out position;V_{y} (fm);#", 500, -50., 50.);
+    hVz = new TH1F("hVz", "z freeze-out position;V_{z} (fm);#", 500, -50., 50.);
+    hPx = new TH1F("hPx", "single particle p_{x};p_{x} (GeV/c);#", 200, -2., 2.);
+    hPy = new TH1F("hPy", "single particle p_{y};p_{y} (GeV/c);#", 200, -2., 2.);
+    hPz = new TH1F("hPz", "single particle p_{z};p_{z} (GeV/c);#", 200, -2., 2.);
+    hP  = new TH1F("hP", "single particle p;p (GeV/c);#", 400, 0., 2.);
+    hPt = new TH1F("hPt", "single particle p_{T};p_{T} (GeV/c);#", 400, 0., 2.);  
+    hEta = new TH1F("hEta", ";#eta;#", 200, -10., 10.);
+    hEtaInt = new TH1F("hEtaInt", ";#eta;#", 200, -10., 10.);  
+    hMSqrVsPt = new TH2F("hMSqrVsPt",
+                         "m^{2} vs p_{T};p_{T}*charge (GeV/c);m^{2} (GeV^{2})",
+                         200, -2., 2., // Pt
+                         100, 0., 2.); // msqr
+    hMomElectron = new TH1F("hMomElectron", "Momentum of e^{#pm};p (GeV/c);#frac{dN}{dp}",
+                            250, 0., 2.5);
+    hMomPion     = new TH1F("hMomPion", "Momentum of #pi^{#pm};p (GeV/c);#frac{dN}{dp}",
+                            250, 0., 2.5);
+    hMomKaon     = new TH1F("hMomKaon", "Momentum of K^{#pm};p (GeV/c);#frac{dN}{dp}",
+                            250, 0., 2.5);
+    hMomProton   = new TH1F("hMomProton", "Momentum of p^{#pm};p (GeV/c);#frac{dN}{dp}",
+                            250, 0., 2.5);
+
+    hRfrVsTfr_pion = new TH2F("hRfrVsTfr_pion", ";R^{#pi^{ch}}_{freeze out} [fm];#tau^{#pi^{ch}}_{freeze out} [fm/c]",
+                              300, 0., 15.,
+                              300, 0., 15.);
+    hRfrVsTfr_kaon = new TH2F("hRfrVsTfr_kaon", ";R^{K^{ch}}_{freeze out} [fm];#tau^{K^{ch}}_{freeze out} [fm/c]",
+                              300, 0., 15.,
+                              300, 0., 15.);
+    hRfrVsTfr_proton = new TH2F("hRfrVsTfr_proton", ";R^{p^{#pm}}_{freeze out} [fm];#tau^{p^{#pm}}_{freeze out} [fm/c]",
+                                300, 0., 15.,
+                                300, 0., 15.);
+
+    hdNdTfr_pion = new TH1F("hdNdTfr_pion", ";#frac{dN}{d#tau^{#pi^{ch}_{freezeout}}};#tau^{#pi^{ch}}_{freezeout} [fm/c]",
+                            1000, 0., 1000.);
+    hdNdTfr_kaon = new TH1F("hdNdTfr_kaon", ";#frac{dN}{d#tau^{K^{ch}_{freezeout}}};#tau^{K^{ch}}_{freezeout} [fm/c]",
+                            1000, 0., 1000.);
+    hdNdTfr_proton = new TH1F("hdNdTfr_proton", ";#frac{dN}{d#tau^{p^{ch}_{freezeout}}};#tau^{p^{ch}}_{freezeout} [fm/c]",
+                              1000, 0., 1000.);
+
+    return kStOk;
 }
 
 
 //_________________
 Int_t StO97DstQAMaker::Finish() {
-  UninitRead();
+    UninitRead();
+
+    mOutFile->Write();
+    mOutFile->Close();
 
-  mOutFile->Write();
-  mOutFile->Close();
-  
-  delete pdgDb;
+    delete pdgDb;
 
-  return kStOk;
+    return kStOk;
 }
 
 
 //_________________
 int StO97DstQAMaker::FillChain(TChain *chain, char *dir,
-			       const char *filter, int maxFiles) {
-  TSystem *gSystem = new TSystem();
-  if (!gSystem) {
-    std::cout << "[StO97DstQAMaker]: can't allocate memory for TSystem" << std::endl;;
-    return 0;
-  }
-  
-  void *gDir = gSystem->OpenDirectory(dir);
-  if (!gDir) {
-    std::cout << "[StO97DstQAMaker]: can't open directory " 
-	      << dir << std::endl;
-    delete gSystem;
-    return 0;
-  }
-  
-  int count = 0;
-  char *file;
-
-  while ((file = (char *)gSystem->GetDirEntry(dir))) {
-    if (strcmp(file, ".") == 0 || strcmp(file, "..") == 0)
-      continue;
-    if (strstr(file, filter)) {
-      char *gFullName = gSystem->ConcatFileName(dir, file);
-      std::cout << "[StO97DstQAMaker]: Adding "
-		<< gFullName << " to the chain" << std::endl;
-      chain->Add(gFullName);
-      delete gFullName;
-      count++;
-      if ((maxFiles > 0) && (count > maxFiles)) break;
+                               const char *filter, int maxFiles) {
+    TSystem *gSystem = new TSystem();
+    if (!gSystem) {
+        std::cout << "[StO97DstQAMaker]: can't allocate memory for TSystem" << std::endl;;
+        return 0;
+    }
+
+    void *gDir = gSystem->OpenDirectory(dir);
+    if (!gDir) {
+        std::cout << "[StO97DstQAMaker]: can't open directory " 
+            << dir << std::endl;
+        delete gSystem;
+        return 0;
+    }
+
+    int count = 0;
+    char *file;
+
+    while ((file = (char *)gSystem->GetDirEntry(dir))) {
+        if (strcmp(file, ".") == 0 || strcmp(file, "..") == 0)
+            continue;
+        if (strstr(file, filter)) {
+            char *gFullName = gSystem->ConcatFileName(dir, file);
+            std::cout << "[StO97DstQAMaker]: Adding "
+                << gFullName << " to the chain" << std::endl;
+            chain->Add(gFullName);
+            delete gFullName;
+            count++;
+            if ((maxFiles > 0) && (count > maxFiles)) break;
+        }
     }
-  }
-  
-  std::cout << "[StO97DstQAMaker]: Added " << count 
-	    << " files to the chain" << std::endl;
-  delete gSystem;
-  return count;
+
+    std::cout << "[StO97DstQAMaker]: Added " << count 
+        << " files to the chain" << std::endl;
+    delete gSystem;
+    return count;
 }
 
 
 //_________________
 int StO97DstQAMaker::FillChain(TChain *chain, const char *fileName,
-			       int maxFiles) {
-  std::ifstream inStream(fileName);
-
-
-  if (!inStream.is_open()) {
-    std::cout << "[StFemtoDstQAMaker]: can't open file list: "
-	      << fileName << std::endl;
-    return 0;
-  }
-  std::cout << "[StFemtoDstQAMaker]: using file list: "
-	    << fileName << std::endl;
-
-  int count = 0;
-  char buf[0xFF];
-  for ( ; inStream.good(); ) {
-    inStream.getline(buf, 0xFF);
-    TString lFileName(buf);
-    if (lFileName.Contains("root")) {
-      chain->Add(buf);
-      count++;
-      if ((maxFiles > 0) && (count > maxFiles)) break;
+                               int maxFiles) {
+    std::ifstream inStream(fileName);
+
+
+    if (!inStream.is_open()) {
+        std::cout << "[StFemtoDstQAMaker]: can't open file list: "
+            << fileName << std::endl;
+        return 0;
+    }
+    std::cout << "[StFemtoDstQAMaker]: using file list: "
+        << fileName << std::endl;
+
+    int count = 0;
+    char buf[0xFF];
+    for ( ; inStream.good(); ) {
+        inStream.getline(buf, 0xFF);
+        TString lFileName(buf);
+        if (lFileName.Contains("root")) {
+            chain->Add(buf);
+            count++;
+            if ((maxFiles > 0) && (count > maxFiles)) break;
+        }
     }
-  }
 
-  std::cout << "[StO97DstQAMaker]: Added " << count
-	    << " files to the chain" << std::endl;
-  return count;
+    std::cout << "[StO97DstQAMaker]: Added " << count
+        << " files to the chain" << std::endl;
+    return count;
 }
 
 
 
 //_________________
 int StO97DstQAMaker::InitRead(string dir, string fileName,
-			      string filter, int maxFiles) {
-  mEventId = 0;
-
-  mChain = new TChain("StO97Dst", "StO97Dst");
-  
-  std::cout << "[StO97DstQAMaker]: Call InitRead" << std::endl;
-  
-  int numFiles = 0;
-  if (!fileName.empty()) {
-    if (strstr(fileName.c_str(), ".lis") ||
-	strstr(fileName.c_str(), ".list")) {
-      numFiles = FillChain(mChain, 
-			   (dir + fileName).c_str(), 
-			   maxFiles);
+                              string filter, int maxFiles) {
+    mEventId = 0;
+
+    mChain = new TChain("StO97Dst", "StO97Dst");
+
+    std::cout << "[StO97DstQAMaker]: Call InitRead" << std::endl;
+
+    int numFiles = 0;
+    if (!fileName.empty()) {
+        if (strstr(fileName.c_str(), ".lis") ||
+            strstr(fileName.c_str(), ".list")) {
+            numFiles = FillChain(mChain, 
+                                 (dir + fileName).c_str(), 
+                                 maxFiles);
+        }
+        else {
+            mChain->Add((dir + fileName).c_str());
+            numFiles++;
+        }
+    } else {
+        numFiles = FillChain(mChain, (char *)dir.c_str(), // not cool
+                             filter.c_str(), maxFiles);
     }
-    else {
-      mChain->Add((dir + fileName).c_str());
-      numFiles++;
-    }
-  } else {
-    numFiles = FillChain(mChain, (char *)dir.c_str(), // not cool
-			 filter.c_str(), maxFiles);
-  }
-  mChain->SetBranchAddress("StO97Event", &mO97Event);
-  mNEvents = (unsigned int)mChain->GetEntries();
-  return numFiles;
+    mChain->SetBranchAddress("StO97Event", &mO97Event);
+    mNEvents = (unsigned int)mChain->GetEntries();
+    return numFiles;
 }
 
 
 //_________________
 void StO97DstQAMaker::UninitRead() {
-  if (mChain)     delete mChain;
-  mO97Event = 0;
-  mChain    = 0;
+    if (mChain)     delete mChain;
+    mO97Event = 0;
+    mChain    = 0;
 }
 
 
 //_________________
 int StO97DstQAMaker::GetNEvents() {
-  if (mChain)
-    return mNEvents;
-  else
-    return -1;
+    if (mChain)
+        return mNEvents;
+    else
+        return -1;
 }
 
 
 //_________________
 Int_t StO97DstQAMaker::Make() {
-  if (!mNEvents) {
-    std::cout << "[StO97DstQAMaker]: no events to read" << std::endl;
-    return 0;
-  }
-
-  mO97Event->Clear();
-  int bytes = mChain->GetEntry(mEventId++);
-
-  if (mNEvents < mEventId) {
-    std::cout << "[StO97DstQAMaker]: EOF" << std::endl;
-    return 0;
-  }
-
-  if (!bytes) {
-    std::cout << "[StO97DstQAMaker]: no event" << std::endl;
-    return 0;
-  }
-
-  if (mO97Event) {
-    int refMultP90 = 0;
-    int refMultP100 = 0;
-    int refMultP120 = 0;
-    int refMultP150 = 0;
-    float imp = mO97Event->GetImpactPar();
-    
-    TClonesArray *tracks = mO97Event->GetTracks();
-    
-    if (tracks) {
-      int nTracks = tracks->GetEntries();
-      for (int i = 0; i < nTracks; i++) {
-	StO97Track *track = (StO97Track *)tracks->At(i);
-	float px = track->GetPx();
-	float py = track->GetPy();
-	float pz = track->GetPz();
-	float pt = sqrt(px*px + py*py);
-	float p = sqrt(px*px + py*py + pz*pz);
-	float massSqr = track->GetMass()*track->GetMass();
-	float vx = track->GetXfr();
-	float vy = track->GetYfr();
-	float vz = track->GetZfr();
-	float eta = 0.5*log((p + pz)/(p - pz));
-	int pdg = track->GetPdgId();
-	
-	TParticlePDG *pdgInfo = pdgDb->GetParticle(pdg);
-
-	hEtaInt->Fill(eta);
-
-	if (pt < 0.075 || fabs(eta) >= 0.5 || pdgInfo->Charge() == 0. || track->GetIsSpec()) continue;
-
-	if (pt >= 0.09)   refMultP90++;
-	if (pt >= 0.10)   refMultP100++;
-	if (pt >= 0.12)   refMultP120++;
-	if (pt >= 0.15)   refMultP150++;
-	
-	if (pt >= 0.085) {
-	  hMSqrVsPt->Fill(pt*pdgInfo->Charge(), massSqr);
-	  hVx->Fill(vx);
-	  hVy->Fill(vy);
-	  hVz->Fill(vz);
-	  hPx->Fill(px);
-	  hPy->Fill(py);
-	  hPz->Fill(pz);
-	  hPt->Fill(pt);
-	  hP->Fill(p);
-	  hEta->Fill(eta);
-	}
-	
-	switch (abs(pdg)) {
-	case 211:
-	case 321:
-	case 2212:
-	  break;
-	} // switch (abs(pdg))
-      } // for (int i = 0; i < nTracks; i++)
-    } // if (tracks)
-    else {
-      std::cout << "[StO97DstQAMaker]: tracks == NULL" << std::endl;
+    if (!mNEvents) {
+        std::cout << "[StO97DstQAMaker]: no events to read" << std::endl;
+        return 0;
+    }
+
+    mO97Event->Clear();
+    int bytes = mChain->GetEntry(mEventId++);
+
+    if (mNEvents < mEventId) {
+        std::cout << "[StO97DstQAMaker]: EOF" << std::endl;
+        return 0;
     }
-    
-    if (refMultP90)
-      hRefMultP90->Fill(refMultP90);
-    if (refMultP100)
-      hRefMultP100->Fill(refMultP100);
-    if (refMultP120)
-      hRefMultP120->Fill(refMultP120);
-    if (refMultP150)
-      hRefMultP150->Fill(refMultP150);
-
-    hImpPar->Fill(imp);
-    //    hImpParVsRefNik->Fill(imp, refMultNik);
-  } // if (mO97Event)
-  
-  return kStOk;
+
+    if (!bytes) {
+        std::cout << "[StO97DstQAMaker]: no event" << std::endl;
+        return 0;
+    }
+
+    if (mO97Event) {
+        int refMultP90 = 0;
+        int refMultP100 = 0;
+        int refMultP120 = 0;
+        int refMultP150 = 0;
+        float imp = mO97Event->GetImpactPar();
+
+        TClonesArray *tracks = mO97Event->GetTracks();
+
+        if (tracks) {
+            int nTracks = tracks->GetEntries();
+            for (int i = 0; i < nTracks; i++) {
+                StO97Track *track = (StO97Track *)tracks->At(i);
+                float px = track->GetPx();
+                float py = track->GetPy();
+                float pz = track->GetPz();
+                float pt = sqrt(px*px + py*py);
+                float p = sqrt(px*px + py*py + pz*pz);
+                float massSqr = track->GetMass()*track->GetMass();
+                float vx = track->GetXfr();
+                float vy = track->GetYfr();
+                float vz = track->GetZfr();
+                float Rfr = TMath::Sqrt(vx*vx + vy*vy);
+                float Tfr = track->GetTfr();
+                float eta = pt > 0.05 /* i.e. pt > 0 */ ? 0.5*log((p + pz)/(p - pz)) : -999.;
+                int   pdg = track->GetPdgId();
+                int   charge = track->GetCharge(pdgDb);
+
+                hEtaInt->Fill(eta);
+
+                if (pt < 0.075 || fabs(eta) >= 0.5 /* STAR TPC acceptance */ ||
+                    charge == 0. /* only charged particles */ ||
+                    track->GetIsSpec()) {
+
+                    continue;
+                }
+
+                if (pt >= 0.09)   refMultP90++;
+                if (pt >= 0.10)   refMultP100++;
+                if (pt >= 0.12)   refMultP120++;
+                if (pt >= 0.15)   refMultP150++;
+
+                if (pt >= 0.1 /* STAR TPC lower momentum limit from STAR TPC NIM */) {
+                    hMSqrVsPt->Fill(pt*charge, massSqr);
+                    hVx->Fill(vx);
+                    hVy->Fill(vy);
+                    hVz->Fill(vz);
+                    hPx->Fill(px);
+                    hPy->Fill(py);
+                    hPz->Fill(pz);
+                    hPt->Fill(pt);
+                    hP->Fill(p);
+                    hEta->Fill(eta);
+                }
+
+                switch (abs(pdg)) {
+                    case 11: // Electron
+                        hMomElectron->Fill(p);
+                        break;
+                    case 211: // Pion
+                        hMomPion->Fill(p);
+                        hdNdTfr_pion->Fill(Tfr);
+                        hRfrVsTfr_pion->Fill(Rfr, Tfr);
+                        break;
+                    case 321: // Kaon
+                        hMomKaon->Fill(p);
+                        hdNdTfr_kaon->Fill(Tfr);
+                        hRfrVsTfr_kaon->Fill(Rfr, Tfr);
+                        break;
+                    case 2212: // Proton
+                        hMomProton->Fill(p);
+                        hdNdTfr_proton->Fill(Tfr);
+                        hRfrVsTfr_proton->Fill(Rfr, Tfr);
+                } // switch (abs(pdg))
+            } // for (int i = 0; i < nTracks; i++)
+        } // if (tracks)
+        else {
+            std::cout << "[StO97DstQAMaker]: tracks == NULL" << std::endl;
+        }
+
+        if (refMultP90)
+            hRefMultP90->Fill(refMultP90);
+        if (refMultP100)
+            hRefMultP100->Fill(refMultP100);
+        if (refMultP120)
+            hRefMultP120->Fill(refMultP120);
+        if (refMultP150)
+            hRefMultP150->Fill(refMultP150);
+
+        hSph->Fill(mO97Event->GetTransverseSphericity());
+        hSph2->Fill(mO97Event->GetTransverseSphericity2());
+
+        hImpPar->Fill(imp);
+				hImpVsRef->Fill(imp, refMultP150);
+    } // if (mO97Event)
+
+    return kStOk;
 }

+ 31 - 12
StRoot/StFemtoDstMaker/StO97DstQAMaker.h

@@ -1,3 +1,10 @@
+/*
+ * READ ME:
+ *   For proper working of this maker you need to link StRoot to your directory.
+ *   For more info please take a look on StO97DstQAMaker.cxx line 37
+ *
+ */
+
 #ifndef STO97DSTQAMAKER_HH
 #define STO97DSTQAMAKER_HH
 
@@ -12,7 +19,6 @@
 #include "StO97Event.h"
 #include <string>
 #include <stdlib.h>
-#include "TParticlePDG.h"
 #include "TDatabasePDG.h"
 
 
@@ -20,11 +26,11 @@
 class StO97DstQAMaker : public StMaker {
  public:
   StO97DstQAMaker(const char *dirName,
-		  const char *fileName,
-		  const char *filter = ".",
-		  int        maxFiles = 1e5);
+                  const char *fileName,
+                  const char *filter = ".",
+                  int        maxFiles = 1e5);
   int GetNEvents();
-  
+
   Int_t Init();
   Int_t Make();
   Int_t Finish();
@@ -53,9 +59,10 @@ class StO97DstQAMaker : public StMaker {
   TH1F *hRefMultP120;
   TH1F *hRefMultP150;
   TH1F *hImpPar;
-  TH2F *hImpParVsRefMNik;
-  TH2F *hImpParVsRefMGreg;
   TH1F *hImp;
+  TH1F *hSph;
+  TH1F *hSph2;
+	TH2F *hImpVsRef;
 
   //
   // Track histograms
@@ -71,15 +78,27 @@ class StO97DstQAMaker : public StMaker {
   TH1F *hPz;
   TH1F *hEta;
   TH1F *hEtaInt;
-  
+  TH1F *hMomElectron;
+  TH1F *hMomPion;
+  TH1F *hMomKaon;
+  TH1F *hMomProton;
+
+  TH2F *hRfrVsTfr_pion;
+  TH2F *hRfrVsTfr_kaon;
+  TH2F *hRfrVsTfr_proton;
+
+  TH1F *hdNdTfr_pion;
+  TH1F *hdNdTfr_kaon;
+  TH1F *hdNdTfr_proton;
+
   int InitRead(string dir, string fileName,
-	       string filter, int maxFiles);
+               string filter, int maxFiles);
   int FillChain(TChain *chain, const char *fileName, 
-		int maxFiles);
+                int maxFiles);
   int FillChain(TChain *chain, char *dir,
-		const char *filter, int maxFiles);
+                const char *filter, int maxFiles);
   void UninitRead();
-  
+
   ClassDef(StO97DstQAMaker, 1)
 };
 

+ 3 - 2
StRoot/StFemtoDstMaker/StO97Event.cxx

@@ -4,9 +4,10 @@
 ClassImp(StO97Event)
 
 //_________________
-StO97Event::StO97Event() : mRunId(0), mEventNumber(0), mNTracks(0),
-  mImpactPar(0.), mEventPlaneRot(0.), mTransverseSphericity(0),
+StO97Event::StO97Event() : mRunId(0), mEventNumber(0),
+  mImpactPar(0.), mEventPlaneRot(0.), mNTracks(0), mTransverseSphericity(0),
   mTransverseSphericity2(0) {
+
   mTracks = new TClonesArray("StO97Track", 1500);
 }
 

+ 0 - 1
StRoot/StFemtoDstMaker/StO97Event.h

@@ -39,7 +39,6 @@ class StO97Event : public TObject {
   //Setters
   void SetRunId(Int_t id)                { (id>0) ? mRunId=(UInt_t)id : mRunId=(UInt_t)(-1 * id); }
   void SetEventNumber(Int_t num)         { (num>0) ? mEventNumber=(UInt_t)num : mEventNumber=(UInt_t)(-1 * num); }
-  void SetNTracks(Int_t nTracks)         { mNTracks = nTracks; }
   void SetImpactPar(Float_t impPar)      { mImpactPar = impPar; }
   void SetEventPlaneRot(Float_t epr)     { mEventPlaneRot = epr; }
   void SetTracks(TClonesArray *tracks)   { mTracks = tracks; }

+ 21 - 26
StRoot/StFemtoDstMaker/StO97Track.cxx

@@ -1,6 +1,4 @@
 #include "StO97Track.h"
-#include <TMath.h>
-#include <iostream>
 
 ClassImp(StO97Track)
 
@@ -21,6 +19,12 @@ StO97Track::StO97Track(const StO97Track &t) {
 }
 
 //________________
+Int_t StO97Track::GetPdgId(const TDatabasePDG *pdgDb)
+{
+  return pdgDb->ConvertIsajetToPdg(mPdgId);
+}
+
+//________________
 StO97Track &StO97Track::operator=(const StO97Track &t) {
 
   //Assignment constructor
@@ -40,6 +44,19 @@ StO97Track &StO97Track::operator=(const StO97Track &t) {
 }
 
 //_________________
+Int_t StO97Track::GetCharge(const TDatabasePDG *pdgDb, bool isajet) {
+  int pdg = isajet ? pdgDb->ConvertIsajetToPdg(mPdgId) : mPdgId;
+  TParticlePDG *pdgInfo = pdgDb->GetParticle(pdg);
+  if (pdgInfo != NULL) {
+    return (Int_t)(pdgInfo->Charge()/3.);
+  }
+  else {
+    std::cout << "[StO97Track WARN]: can't get TPartcilePDG, using mPdgId > 0 ? 1 : -1\n";
+    return GetCharge();
+  }
+}
+
+//_________________
 StO97Track::~StO97Track() {
   /* empty */
 }
@@ -57,23 +74,10 @@ void StO97Track::SetId(Int_t id) {
   }
 }
 
-//_________________
-void StO97Track::SetPdgId(Int_t id) {
-
-  //Saving PdgId of the track
-  if( ( id < std::numeric_limits<unsigned short>::min() ) ||
-      ( id > std::numeric_limits<unsigned short>::max() ) ) {
-    mPdgId = 0;
-  }
-  else {
-    mPdgId = id;
-  }
-}
-
 //________________
 void StO97Track::SetTrack(Int_t id, Int_t pdgId, 
-			  Float_t px, Float_t py, Float_t pz, Float_t mass,
-			  Float_t xFr, Float_t yFr, Float_t zFr, Float_t tFr) {
+                          Float_t px, Float_t py, Float_t pz, Float_t mass,
+                          Float_t xFr, Float_t yFr, Float_t zFr, Float_t tFr) {
   
   //Set all track information
   SetId(id);
@@ -136,12 +140,3 @@ Float_t StO97Track::GetRapidity() {
 Float_t StO97Track::GetY() {
   return GetRapidity();
 }
-
-//_________________
-void StO97Track::Print() {
-  std::cout << "Track params: " <<
-	    << "id: " << mId << " pdgId: " << mPdgId
-	    << " px: " << mPx << " py: " << mPy << " pz: " << mPz
-	    << " mass: " << mMass << " xFr: " << mXfr << " yFr: " << mYfr
-	    << " zFr: " << mZfr << "tFr: " << mTfr << std::endl;
-}

+ 17 - 3
StRoot/StFemtoDstMaker/StO97Track.h

@@ -3,11 +3,17 @@
 //  when any changes in the codes are done
 //  Grigory Nigmatkulov: 2016/12/15
 //
+
 #ifndef STOSCAR97TRACK_H
 #define STOSCAR97TRACK_H
 
 #include <TObject.h>
 #include <limits>
+#include <iostream>
+#include <TMath.h>
+#include "TParticlePDG.h"
+#include "TDatabasePDG.h"
+
 
 //_________________
 class StO97Track : public TObject {
@@ -26,7 +32,9 @@ class StO97Track : public TObject {
   //Getters
   Int_t   GetId() const     { return (Int_t)mId; }
   Int_t   GetPdgId() const  { return mPdgId; }
+  Int_t   GetPdgId(const TDatabasePDG *pdgDb);  
   Int_t   GetCharge() const { return (mPdgId>0) ? 1 : -1; }
+  Int_t   GetCharge(const TDatabasePDG *pdgDb, bool isajet = false);
   Float_t GetPx() const     { return mPx; }
   Float_t GetPy() const     { return mPy; }
   Float_t GetPz() const     { return mPz; }
@@ -48,7 +56,7 @@ class StO97Track : public TObject {
 
   //Setters
   void SetId(Int_t id);
-  void SetPdgId(Int_t pdgId);
+  void SetPdgId(Int_t pdgId) { mPdgId = pdgId; }
   void SetPx(Float_t px)     { mPx = px; }
   void SetPy(Float_t py)     { mPy = py; }
   void SetPz(Float_t pz)     { mPz = pz;}
@@ -63,8 +71,14 @@ class StO97Track : public TObject {
 		Float_t xFr, Float_t yFr, Float_t zFr, Float_t tFr);
 
   //Print track parameters
-  void Print();
-  
+  virtual void Print(Option_t *option="") const {
+      std::cout << "Track params: "
+          << "id: " << mId << " pdgId: " << mPdgId
+          << " px: " << mPx << " py: " << mPy << " pz: " << mPz
+          << " mass: " << mMass << " xFr: " << mXfr << " yFr: " << mYfr
+          << " zFr: " << mZfr << "tFr: " << mTfr << std::endl;
+  }
+
  private:
 
   UShort_t mId;     // Number of the particle

File diff suppressed because it is too large
+ 3797 - 0
StRoot/StFemtoDstMaker/pdg_table.txt


+ 592 - 0
StRoot/StFemtoDstMakerNORP/StFemtoDstInclusiveSelector.cxx

@@ -0,0 +1,592 @@
+#include "StFemtoDstInclusiveSelector.h"
+#include "StBTofHeader.h"
+#include "TBranch.h"
+#include "StPhysicalHelixD.hh"
+
+ClassImp(StFemtoDstInclusiveSelector)
+
+// Set maximum file size to 1.9 GB (Root has a 2GB limit)
+#define MAXFILESIZE 1900000000
+
+//_________________
+StFemtoDstInclusiveSelector::StFemtoDstInclusiveSelector(StMuDstMaker *muMaker,
+                                                         const Char_t *oFileName) {
+
+  std::cout << "StFemtoDstInclusiveSelector::Constructor - Creating an instance...";
+  mOutFileName = oFileName;
+  mMuDstMaker = muMaker;
+  mMuDst = mMuDstMaker->muDst();
+  mRefMultCorrUtil = NULL;
+  mCompression = 9;
+
+  mCollisionTypeAuAu = false; // pp: false; AuAu: true
+  mLoopAllVert = false;       // Not first primary vertex does not contain tracks with fast detectors (TOF)
+
+  mEventIsGood = false;
+  mNEventsIn = 0;
+  mNEventsPassed = 0;
+
+  mAuAuCorrectZdc = true;
+
+  mIsGoodTrack = false;
+  mCurrentTrigger = 0;
+
+  //Initialize event cut variables
+  mPrimVtxZ[0] = -70.;
+  mPrimVtxZ[1] = 70.;
+  mPrimVtxR[0] = -1.;
+  mPrimVtxR[1] = 10.;
+  mPrimVtxVpdVzDiff[0] = -5.;
+  mPrimVtxVpdVzDiff[1]= 5.;
+  mPrimVtxXShift = 0.;
+  mPrimVtxYShift = 0.;
+
+  //Initialize single-particle cut variables
+  mTrackMomentum[0] = 0.1;
+  mTrackMomentum[1] = 1.7;
+  mTrackDca[0] = 0.;
+  mTrackDca[1] = 5.;
+  mTrackDcaGlobal[0] = 0.;
+  mTrackDcaGlobal[1] = 5.;
+  mTrackNHits[0] = 14;
+  mTrackNHits[1] = 50;
+
+  mTrackEta[0] = -1.;
+  mTrackEta[1] = 1.;
+  mTrackFlag[0] = 0;
+  mTrackFlag[1] = 1000;
+
+  mSelectPions = true;
+  mSelectTpcPions = true;
+  mPionTpcMomentum[0] = 0.15;      mPionTpcMomentum[1] = 0.65;
+  mPionTpcNSigmaElectron[0]= -2.; mPionTpcNSigmaElectron[1]= 2.;
+  mPionTpcNSigmaPion[0]= -3.;      mPionTpcNSigmaPion[1]= 3.;
+  mPionTpcNSigmaKaon[0]= -1.5;     mPionTpcNSigmaKaon[1]= 1.5;
+  mPionTpcNSigmaProton[0]= -3.;   mPionTpcNSigmaProton[1]= 3.;
+  mSelectTofPions = true;
+  mPionTofMomentum[0] = 0.15;      mPionTofMomentum[1] = 1.75;
+  mPionMassSqr[0] = -0.01;         mPionMassSqr[1] = 0.09;
+
+  mSelectKaons = true;
+  mSelectTpcKaons = true;
+  mKaonTpcMomentum[0] = 0.15;      mKaonTpcMomentum[1] = 0.65;
+  mKaonTpcNSigmaElectron[0]= -2.; mKaonTpcNSigmaElectron[1]= 2.;
+  mKaonTpcNSigmaPion[0]= -1.5;     mKaonTpcNSigmaPion[1]= 1.5;
+  mKaonTpcNSigmaKaon[0]= -3.;      mKaonTpcNSigmaKaon[1]= 3.;
+  mKaonTpcNSigmaProton[0]= -3.;   mKaonTpcNSigmaProton[1]= 3.;
+  mSelectTofKaons = true;
+  mKaonTofMomentum[0] = 0.15;      mKaonTofMomentum[1] = 1.75;
+  mKaonMassSqr[0] = 0.16;          mKaonMassSqr[1] = 0.35;
+
+  mSelectProtons = false;
+  mSelectTpcProtons = true;
+  mProtonTpcMomentum[0] = 0.15;      mProtonTpcMomentum[1] = 1.2;
+  mProtonTpcNSigmaElectron[0]= -3.; mProtonTpcNSigmaElectron[1]= 3.;
+  mProtonTpcNSigmaPion[0]= -3.;     mProtonTpcNSigmaPion[1]= 3.;
+  mProtonTpcNSigmaKaon[0]= -3.;     mProtonTpcNSigmaKaon[1]= 3.;
+  mProtonTpcNSigmaProton[0]= -3.;    mProtonTpcNSigmaProton[1]= 3.;
+  mSelectTofProtons = true;
+  mProtonTofMomentum[0] = 0.15; mProtonTofMomentum[1] = 1.75;
+  mProtonMassSqr[0] = 0.65; mProtonMassSqr[1] = 1.05;
+
+  std::cout << "\t[DONE]" << std::endl;
+}
+
+//_________________
+StFemtoDstInclusiveSelector::~StFemtoDstInclusiveSelector() {
+  /* do nothing */
+}
+
+//_________________
+Int_t StFemtoDstInclusiveSelector::Init() {
+
+  std::cout << "StFemtoDstInclusiveSelector::Init - Initializing the maker"
+      << std::endl 
+      << "Creating output file: " << mOutFileName;
+
+  mFemtoEvent = new StFemtoEvent();
+  mOutFile = new TFile(mOutFileName, "RECREATE");
+  mOutFile->SetCompressionLevel(mCompression);
+  std::cout << "\t[DONE]" << std::endl;
+
+  mTree = new TTree("StFemtoDst","StFemtoDst");
+  mTree->SetMaxTreeSize(MAXFILESIZE);
+  //mTree->SetAutoSave(1000000);
+  mTree->Branch("StFemtoEvent","StFemtoEvent", &mFemtoEvent);
+
+  if(mCollisionTypeAuAu) { //For AuAu only
+    if(!mRefMultCorrUtil) {
+      std::cout << "StFemtoDstInclusiveSelector::Init - Initializing StRefMultCorr...";
+      mRefMultCorrUtil = new StRefMultCorr("refmult");
+      std::cout << "\t[DONE]" << endl;
+    }
+  }
+
+  mNBytes = 0;
+  std::cout << "StFemtoDstInclusiveSelector::Init - Initialization has been finished"
+      << std::endl;
+  return StMaker::Init();
+}
+
+//_________________
+void StFemtoDstInclusiveSelector::Clear(Option_t *option) {
+  StMaker::Clear(option);
+}
+
+//_________________
+Int_t StFemtoDstInclusiveSelector::Make() {
+
+  // Get MuEvent from the MuDst
+  mNEventsIn++;
+  StMuEvent *mMuEvent = mMuDst->event();
+  if (!AcceptTrigger(mMuEvent))  { // If trigger the trigger was not found
+    return kStOk;
+  }
+  if (mMuEvent->refMult() < 0) { // Multiplicity can't be negative
+    return kStOk;
+  }
+
+  UShort_t mNVertices = mLoopAllVert ? mMuDst->numberOfPrimaryVertices() : 1;
+  Short_t mNPrimTracks = mMuDst->numberOfPrimaryTracks();
+  Short_t mNGlobTracks = mMuDst->numberOfGlobalTracks();
+
+  //Reasonable amount of tracks
+  if (mNPrimTracks < 0 || mNPrimTracks > 10000) {
+    return kStOk;
+  }
+
+  Bool_t mEventIsGood = false;
+  Float_t mVpdVz = 0.;
+  StBTofHeader *mTofHeader = mMuDst->btofHeader();
+
+  // Centrality variables
+  Int_t mCent16;
+  Double_t mCentWeight, mRefMultCorrVal;
+
+  // If pp collisions
+  if (!mCollisionTypeAuAu) {
+    mCent16 = -1;
+    mCentWeight = -999.;
+  }
+
+  Float_t mRanking = 0;
+  StMuPrimaryVertex *mPrimVertex = NULL;
+  StMuTrack *mPrimTrack = NULL;
+  StMuTrack *mGlobTrack = NULL;
+  Float_t mMassSqr = -999.;
+  Float_t mBeta = -999.;
+  Bool_t mIsTofTrack = false;
+  Bool_t mIsTpcPion = false;
+  Bool_t mIsTofPion = false;
+  Bool_t mIsTpcKaon = false;
+  Bool_t mIsTofKaon = false;
+  Bool_t mIsTpcProton = false;
+  Bool_t mIsTofProton = false;
+
+  //Calculate the amount of the pile-up vertices
+  UShort_t mNVtxInBunchCross = 0;
+  for(UShort_t iVert=0; iVert<mMuDst->numberOfPrimaryVertices(); iVert++) {
+
+    mPrimVertex = mMuDst->primaryVertex(iVert);
+    if(!mPrimVertex || !mTofHeader) continue;
+    if(mPrimVertex->ranking() <= 0.) continue;
+
+    if(TMath::Abs(mPrimVertex->position().x()) < 1e-5 &&
+       TMath::Abs(mPrimVertex->position().y()) < 1e-5 &&
+       TMath::Abs(mPrimVertex->position().z()) < 1e-5) {
+      continue;
+    }
+    if(TMath::Abs( mPrimVertex->position().z() - mTofHeader->vpdVz()) < 5. ) {
+      mNVtxInBunchCross++;
+    }
+  } //for(unsigned short iVert=0; iVert<mNVertices; iVert++)
+
+  //There should be only one primary vertex that caused the trigger
+  if(mNVtxInBunchCross > 1) {
+    return kStOk;
+  }
+
+  // Vertex loop
+  for (UShort_t iVert = 0; iVert < mNVertices; iVert++) {
+
+    mEventIsGood = false;
+    mPrimVertex = mMuDst->primaryVertex(iVert);
+    if (!mPrimVertex) {
+      //std::cout << "No primary vertices in the event. Skip..." << std::endl;
+      continue;
+    }
+
+    mRanking = mPrimVertex->ranking();
+    //Positive ranking only
+    if (mPrimVertex->ranking() <= 0) continue;
+
+    //Not (0,0,0) position of the primary vertex
+    if (TMath::Abs(mPrimVertex->position().x()) < 1e-5 &&
+        TMath::Abs(mPrimVertex->position().y()) < 1e-5 &&
+        TMath::Abs(mPrimVertex->position().z()) < 1e-5 ) continue;
+
+    StThreeVectorF mVertPosition = mPrimVertex->position();    
+    if(mTofHeader) { //In case of the absence of the vpd information
+      mVpdVz = mTofHeader->vpdVz();
+    }
+    else {
+      mVpdVz = mVertPosition.z();
+    }
+
+    // Event cut
+    if (!AcceptPrimVtx(mVertPosition, mVpdVz)) continue;
+    mEventIsGood = true;
+
+    // Set StRefMultCorr
+    if (mCollisionTypeAuAu) { // Only for AuAu collisions
+      mRefMultCorrUtil->init(mMuEvent->runId());
+      //Check for the bad run
+      if (mRefMultCorrUtil->isBadRun(mMuEvent->runId())) {
+        std::cout << "StRefMultCorrUtil::IsBadRun" << std::endl;
+        continue;
+      }
+      if (mAuAuCorrectZdc == true) {
+        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
+                                    mVertPosition.z(),
+                                    mMuEvent->runInfo().zdcCoincidenceRate());
+      }
+      else {
+        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
+                                    mVertPosition.z());
+      }
+      mCent16 = mRefMultCorrUtil->getCentralityBin16();
+      if (mCent16 < 0) continue;   //Some AuAu collisions have -1 value
+      mCentWeight = mRefMultCorrUtil->getWeight();
+      mRefMultCorrVal = mRefMultCorrUtil->getRefMultCorr();
+    }
+
+    // Write event information to the FemtoDst
+    mNEventsPassed++;
+    mFemtoEvent->SetEventId( mMuEvent->eventId() );
+    mFemtoEvent->SetRunId( mMuEvent->runId() );
+    mFemtoEvent->SetCollisionType( mCollisionTypeAuAu );
+    mFemtoEvent->SetRefMult( mMuEvent->refMult() );
+    mFemtoEvent->SetRefMultCorr( mRefMultCorrVal );
+    mFemtoEvent->SetRefMultCorrWeight( mCentWeight );