Nikita Ermakov 2 years ago
parent
commit
c272a9f268
99 changed files with 13681 additions and 875 deletions
  1. 119 107
      StRoot/StFemtoDstMaker/StFemtoDstInclusiveSelector.cxx
  2. 1 0
      StRoot/StFemtoDstMaker/StFemtoDstInclusiveSelector.h
  3. 172 162
      StRoot/StFemtoDstMaker/StFemtoDstMaker.cxx
  4. 1 0
      StRoot/StFemtoDstMaker/StFemtoDstMaker.h
  5. 3 4
      StRoot/StFemtoDstMaker/StFemtoDstPythia6Maker.cxx
  6. 207 167
      StRoot/StFemtoDstMaker/StFemtoDstQAMaker.cxx
  7. 19 7
      StRoot/StFemtoDstMaker/StFemtoDstQAMaker.h
  8. 11 0
      StRoot/StFemtoDstMaker/StFemtoEvent.cxx
  9. 8 4
      StRoot/StFemtoDstMaker/StFemtoEvent.h
  10. 1 1
      StRoot/StFemtoDstMaker/StFemtoTrack.cxx
  11. 3 3
      StRoot/StFemtoDstMaker/StFemtoTrack.h
  12. 106 86
      StRoot/StFemtoDstMaker/StHbtFemtoDstReader.cxx
  13. 3 2
      StRoot/StFemtoDstMaker/StHbtFemtoDstReader.h
  14. 8 8
      StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_MomCons.cxx
  15. 1 1
      StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_MomCons.h
  16. 1 0
      StRoot/StHbtMaker/Infrastructure/StHbtManager.cxx
  17. 848 0
      StRoot/StMuDstQAMaker/StMuDstMinvMaker.cxx
  18. 281 0
      StRoot/StMuDstQAMaker/StMuDstMinvMaker.h
  19. 9 7
      StRoot/StMuDstQAMaker/StMuDstQAMaker.cxx
  20. BIN
      StRoot/StPicoDstMaker/Diagram1.dia.autosave
  21. BIN
      StRoot/StPicoDstMaker/Diagram1.dia.png
  22. 35 0
      StRoot/StPicoDstMaker/StPicoArrays.cxx
  23. 33 0
      StRoot/StPicoDstMaker/StPicoArrays.h
  24. 40 0
      StRoot/StPicoDstMaker/StPicoBTOWHit.cxx
  25. 32 0
      StRoot/StPicoDstMaker/StPicoBTOWHit.h
  26. 35 0
      StRoot/StPicoDstMaker/StPicoBTofHit.cxx
  27. 31 0
      StRoot/StPicoDstMaker/StPicoBTofHit.h
  28. 65 0
      StRoot/StPicoDstMaker/StPicoBTofPidTraits.cxx
  29. 38 0
      StRoot/StPicoDstMaker/StPicoBTofPidTraits.h
  30. 175 0
      StRoot/StPicoDstMaker/StPicoConstants.cxx
  31. 114 0
      StRoot/StPicoDstMaker/StPicoConstants.h
  32. 247 0
      StRoot/StPicoDstMaker/StPicoCut.cxx
  33. 28 0
      StRoot/StPicoDstMaker/StPicoCut.h
  34. 146 0
      StRoot/StPicoDstMaker/StPicoDst.cxx
  35. 81 0
      StRoot/StPicoDstMaker/StPicoDst.h
  36. 243 0
      StRoot/StPicoDstMaker/StPicoDstHbtReader.cxx
  37. 74 0
      StRoot/StPicoDstMaker/StPicoDstHbtReader.h
  38. 1008 0
      StRoot/StPicoDstMaker/StPicoDstMaker.cxx
  39. BIN
      StRoot/StPicoDstMaker/StPicoDstMaker.cxx.flc
  40. 178 0
      StRoot/StPicoDstMaker/StPicoDstMaker.h
  41. 298 0
      StRoot/StPicoDstMaker/StPicoDstQA.cxx
  42. 115 0
      StRoot/StPicoDstMaker/StPicoDstQA.h
  43. 60 0
      StRoot/StPicoDstMaker/StPicoEmcPidTraits.cxx
  44. 57 0
      StRoot/StPicoDstMaker/StPicoEmcPidTraits.h
  45. 41 0
      StRoot/StPicoDstMaker/StPicoEmcTrigger.cxx
  46. 32 0
      StRoot/StPicoDstMaker/StPicoEmcTrigger.h
  47. 319 0
      StRoot/StPicoDstMaker/StPicoEvent.cxx
  48. 206 0
      StRoot/StPicoDstMaker/StPicoEvent.h
  49. 56 0
      StRoot/StPicoDstMaker/StPicoMtdHit.cxx
  50. 40 0
      StRoot/StPicoDstMaker/StPicoMtdHit.h
  51. 54 0
      StRoot/StPicoDstMaker/StPicoMtdPidTraits.cxx
  52. 47 0
      StRoot/StPicoDstMaker/StPicoMtdPidTraits.h
  53. 103 0
      StRoot/StPicoDstMaker/StPicoMtdTrigger.cxx
  54. 35 0
      StRoot/StPicoDstMaker/StPicoMtdTrigger.h
  55. 208 0
      StRoot/StPicoDstMaker/StPicoTrack.cxx
  56. 150 0
      StRoot/StPicoDstMaker/StPicoTrack.h
  57. 41 0
      StRoot/StPicoDstMaker/StPicoTrigger.cxx
  58. 32 0
      StRoot/StPicoDstMaker/StPicoTrigger.h
  59. 93 0
      StRoot/StPicoDstMaker/StPicoUtilities.h
  60. 150 0
      StRoot/StPicoDstMaker/StPicoV0.cxx
  61. 47 0
      StRoot/StPicoDstMaker/StPicoV0.h
  62. 12 0
      StRoot/StPicoDstMaker/game.000
  63. 3 0
      StRoot/StPicoDstMaker/log.000
  64. 1 0
      StRoot/StPicoDstMaker/players.dat
  65. BIN
      StRoot/StPicoDstMaker/sum_3_0.root
  66. 158 0
      TriggerListsShort.txt
  67. 182 125
      macros/HBT/expKaonHBT_pp200_FemtoDst_1D.C
  68. 501 0
      macros/HBT/expKaonHBT_pp200_FemtoDst_1D_v0.C
  69. 527 0
      macros/HBT/expKaonHBT_pp200_PicoDst_1D.C
  70. 169 119
      macros/HBT/expKaonHBT_pp500_FemtoDst_1D.C
  71. 527 0
      macros/HBT/expKaonHBT_pp500_PicoDst_1D.C
  72. 514 0
      macros/HBT/expKaonHBT_pp500_PicoDst_1D_test.C
  73. 107 0
      macros/Mu2Pico.C
  74. 126 0
      macros/MuDstQA/pau2015QAnia.C
  75. 86 0
      macros/PicoDstQA.C
  76. 15 15
      macros/femtoDst/FemtoDstMaker_kaons_pp200_y2015.C
  77. 220 0
      macros/femtoDst/FemtoDstMaker_mesons_pp200_y2015.C
  78. 37 22
      macros/femtoDst/FemtoDstMaker_mesons_pp500_y2011.C
  79. 104 0
      macros/femtoDst/FemtoDstQA_pp200_y2015.C
  80. 31 5
      macros/femtoDst/FemtoDstQA_pp500_y2011.C
  81. 92 0
      macros/picoDst/pp200_run12_QA.C
  82. 92 0
      macros/picoDst/pp500_run11_QA.C
  83. 2096 0
      processing/pp/DataProcessing.cxx
  84. 370 0
      processing/pp/DataProcessing.h
  85. 137 0
      processing/pp/FemtoUtils.h
  86. 596 0
      processing/pp/HbtFitter.C
  87. 197 0
      processing/pp/HbtFitter.h
  88. 33 0
      processing/pp/Processor.C
  89. 2 0
      sbin/example.list
  90. 28 0
      sbin/genJobs.sh
  91. 23 0
      sbin/hpssMakeList.pl
  92. 10 8
      sbin/listMaker.csh
  93. 2 0
      sbin/out.list
  94. 13 5
      xml/HBT/expKaonHBT_pp200_FemtoDst_1D.xml
  95. 14 6
      xml/HBT/expKaonHBT_pp500_FemtoDst_1D.xml
  96. 28 0
      xml/femtoDst/FemtoDstMaker_mesons_pp200_y2015.xml
  97. 12 9
      xml/femtoDst/FemtoDstMaker_mesons_pp500_y2011.xml
  98. 23 0
      xml/femtoDst/FemtoDstQA_pp200_y2015.xml
  99. 4 2
      xml/femtoDst/FemtoDstQA_pp500_y2011.xml

+ 119 - 107
StRoot/StFemtoDstMaker/StFemtoDstInclusiveSelector.cxx

@@ -277,7 +277,7 @@ Int_t StFemtoDstInclusiveSelector::Make() {
     mFemtoEvent->SetCent16( mCent16 );
     mFemtoEvent->SetRefMultPos( mMuEvent->refMultPos() );
     mFemtoEvent->SetRefMultNeg( mMuEvent->refMultNeg() );
-    mFemtoEvent->SetBTofTrayMultiplicity( mMuEvent->btofTrayMultiplicity() );
+    mFemtoEvent->SetNumberOfBTofHit( mMuDst->numberOfBTofHit() ); //mMuEvent->btofTrayMultiplicity() );
     //mFemtoEvent->SetZDCe( mMuEvent->zdcTriggerDetector().adc(4) );
     //mFemtoEvent->SetZDCw( mMuEvent->zdcTriggerDetector().adc(0) );
     mFemtoEvent->SetNumberOfPrimaryTracks( mNPrimTracks );
@@ -290,6 +290,7 @@ Int_t StFemtoDstInclusiveSelector::Make() {
     mFemtoEvent->SetTriggerId(mCurrentTrigger);
     mFemtoEvent->SetPrimaryVertexRanking(mRanking);
 
+    unsigned int refMult2 = 0;
     //Primary track loop
     for(UShort_t iTrk = 0; iTrk < mNPrimTracks; iTrk++) {
 
@@ -297,6 +298,7 @@ Int_t StFemtoDstInclusiveSelector::Make() {
       mGlobTrack = mMuDst->globalTracks(mPrimTrack->index2Global());
 
       //Track cut
+      if ( AcceptRefMult2(mPrimTrack) ) refMult2++;
       if ( !AcceptTrack(mPrimTrack, iVert) ) continue;
 
       mIsTofTrack = false;
@@ -307,78 +309,79 @@ Int_t StFemtoDstInclusiveSelector::Make() {
       mIsTofKaon = false;
       mIsTpcProton = false;
       mIsTofProton = false;
-      
+
       mBeta = -999.;
       mMassSqr = -999.;
 
       // Check if track has ToF signal
       if(mIsTofTrack) {
-	mBeta = mGlobTrack->btofPidTraits().beta();
-	mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
+        mBeta = mGlobTrack->btofPidTraits().beta();
+        mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
       }
 
       // Pion selection
       if (mSelectPions) {
-	if (mSelectTpcPions) 
-	  mIsTpcPion = IsTpcPion(mPrimTrack);
-	if (mSelectTofPions && mIsTofTrack)
-	  mIsTofPion = IsTofPion(mMassSqr, mGlobTrack->p().mag());
+        if (mSelectTpcPions) 
+          mIsTpcPion = IsTpcPion(mPrimTrack);
+        if (mSelectTofPions && mIsTofTrack)
+          mIsTofPion = IsTofPion(mMassSqr, mGlobTrack->p().mag());
       }
-      
+
       // Kaon selection
       if (mSelectKaons) {
-	if (mSelectTpcKaons)
-	  mIsTpcKaon = IsTpcKaon(mPrimTrack);
-	if (mSelectTofKaons && mIsTofTrack)
-	  mIsTofKaon = IsTofKaon(mMassSqr, mGlobTrack->p().mag());
+        if (mSelectTpcKaons)
+          mIsTpcKaon = IsTpcKaon(mPrimTrack);
+        if (mSelectTofKaons && mIsTofTrack)
+          mIsTofKaon = IsTofKaon(mMassSqr, mGlobTrack->p().mag());
       }
-      
+
       // Proton selection
       if (mSelectProtons) {
-	if (mSelectTpcProtons)
-	  mIsTpcProton = IsTpcProton(mPrimTrack);
-	if (mSelectTofProtons && mIsTofTrack)
-	  mIsTofProton = IsTofProton(mMassSqr, mGlobTrack->p().mag());
+        if (mSelectTpcProtons)
+          mIsTpcProton = IsTpcProton(mPrimTrack);
+        if (mSelectTofProtons && mIsTofTrack)
+          mIsTofProton = IsTofProton(mMassSqr, mGlobTrack->p().mag());
       }
 
       // Write StFemtoTrack
       if (mIsTpcPion || mIsTofPion || mIsTpcKaon || mIsTofKaon || 
-	  mIsTpcProton || mIsTofProton) {
-
-	mFTrack = mFemtoEvent->AddFemtoTrack();
-
-	mFTrack->SetId( mPrimTrack->id() );
-	mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
-	mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss(kTpcId) );
-	mFTrack->SetNHitsFit( mPrimTrack->nHitsFit(kTpcId) );
-	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->SetChi2( mPrimTrack->chi2() );
-	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. );
-	}
+          mIsTpcProton || mIsTofProton) {
+
+        mFTrack = mFemtoEvent->AddFemtoTrack();
+
+        mFTrack->SetId( mPrimTrack->id() );
+        mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
+        mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss(kTpcId) );
+        mFTrack->SetNHitsFit( mPrimTrack->nHitsFit(kTpcId) );
+        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->SetChi2( mPrimTrack->chi2() );
+        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++)
-    
+
+    mFemtoEvent->SetRefMult2(refMult2);
     //
     // If event is good then write it
     //
@@ -392,6 +395,15 @@ Int_t StFemtoDstInclusiveSelector::Make() {
 }
 
 //_________________
+Bool_t StFemtoDstInclusiveSelector::AcceptRefMult2(StMuTrack *trk) {
+  return (trk->flag() > 0 &&
+          trk->charge() != 0 &&
+          trk->nHitsFit() >= 10 &&
+          fabs(trk->eta()) < 1.0 &&
+          trk->dca().mag() < 3.0);
+}
+
+//_________________
 Bool_t StFemtoDstInclusiveSelector::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
@@ -404,9 +416,9 @@ Bool_t StFemtoDstInclusiveSelector::AcceptTrigger(StMuEvent *muEvent) {
   else {
     for (UInt_t iTrg = 0; iTrg < mTriggerIdCollection.size(); iTrg++) {
       if(muEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) {
-	mIsGoodTrigger = true;
-	mCurrentTrigger |= (1<<iTrg);
-	break;
+        mIsGoodTrigger = true;
+        mCurrentTrigger |= (1<<iTrg);
+        break;
       }
     } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
   }
@@ -416,7 +428,7 @@ Bool_t StFemtoDstInclusiveSelector::AcceptTrigger(StMuEvent *muEvent) {
 
 //_________________
 Bool_t StFemtoDstInclusiveSelector::AcceptPrimVtx(StThreeVectorF vtxPos, 
-						  Float_t vpdVz) {
+                                                  Float_t vpdVz) {
   Bool_t mIsGoodVtx = false;
   Float_t mVtxX = vtxPos.x() - mPrimVtxXShift;
   Float_t mVtxY = vtxPos.y() - mPrimVtxYShift;
@@ -424,14 +436,14 @@ Bool_t StFemtoDstInclusiveSelector::AcceptPrimVtx(StThreeVectorF vtxPos,
   Float_t vpdDiff = vtxPos.z() - vpdVz;
 
   if( vtxPos.z() >= mPrimVtxZ[0] &&
-      vtxPos.z() <= mPrimVtxZ[1] &&
-      vtxR >= mPrimVtxR[0] &&
-      vtxR <= mPrimVtxR[1] &&
-      vpdDiff >= mPrimVtxVpdVzDiff[0] &&
-      vpdDiff <= mPrimVtxVpdVzDiff[1]) {
+     vtxPos.z() <= mPrimVtxZ[1] &&
+     vtxR >= mPrimVtxR[0] &&
+     vtxR <= mPrimVtxR[1] &&
+     vpdDiff >= mPrimVtxVpdVzDiff[0] &&
+     vpdDiff <= mPrimVtxVpdVzDiff[1]) {
     mIsGoodVtx = true;
   }
-    
+
   return mIsGoodVtx;
 }
 
@@ -440,24 +452,24 @@ Bool_t StFemtoDstInclusiveSelector::AcceptTrack(StMuTrack *trk, UShort_t vtxInd)
 
   Bool_t mIsGoodTrack = false;
   if( trk->vertexIndex() == vtxInd &&
-      trk->momentum().mag() >= mTrackMomentum[0] && trk->momentum().mag() <= mTrackMomentum[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]) {
+     trk->momentum().mag() >= mTrackMomentum[0] && trk->momentum().mag() <= mTrackMomentum[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]) {
 
     mIsGoodTrack = true;
   }
-      
+
   return mIsGoodTrack;
 }
 
 //_________________
 Bool_t StFemtoDstInclusiveSelector::IsTofTrack(StMuTrack *trk) { //Only for globals
-  
+
   Bool_t mIsTofHit = false;
   if(trk->btofPidTraits().beta() > 0 &&
      trk->btofPidTraits().timeOfFlight() > 0) {
@@ -475,12 +487,12 @@ Int_t StFemtoDstInclusiveSelector::Finish() {
   mOutFile->Close();
   delete mTree;
   delete mOutFile;
-  
+
   std::cout << "*********************************************" << std::endl
-	    << "StFemtoDstInclusiveSelector has been finished" << std::endl
-	    << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
-	    << "\t nEventsProcessed: " << mNEventsIn << std::endl
-	    << "**********************************************" << std::endl;
+      << "StFemtoDstInclusiveSelector has been finished" << std::endl
+      << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
+      << "\t nEventsProcessed: " << mNEventsIn << std::endl
+      << "**********************************************" << std::endl;
 
   return kStOk;
 }
@@ -492,16 +504,16 @@ Bool_t StFemtoDstInclusiveSelector::IsTpcPion(StMuTrack *pTrk) {
   Float_t nSPion = pTrk->nSigmaPion();
   Float_t nSKaon = pTrk->nSigmaKaon();
   Float_t nSProton = pTrk->nSigmaProton();
-  
+
   return ( (p >= mPionTpcMomentum[0]) && (p <= mPionTpcMomentum[1]) &&
-	   ( (nSPion >= mPionTpcNSigmaPion[0]) &&
-	     (nSPion <= mPionTpcNSigmaPion[1]) ) &&
-	   ( (nSElectron <= mPionTpcNSigmaElectron[0]) ||
-	     (nSElectron >= mPionTpcNSigmaElectron[1]))  &&
-	   ( (nSKaon <= mPionTpcNSigmaKaon[0]) ||
-	     (nSKaon >= mPionTpcNSigmaKaon[1]) ) &&
-	   ( (nSProton <= mPionTpcNSigmaProton[0]) ||
-	     (nSProton >= mPionTpcNSigmaProton[1]) ) );
+          ( (nSPion >= mPionTpcNSigmaPion[0]) &&
+           (nSPion <= mPionTpcNSigmaPion[1]) ) &&
+          ( (nSElectron <= mPionTpcNSigmaElectron[0]) ||
+           (nSElectron >= mPionTpcNSigmaElectron[1]))  &&
+          ( (nSKaon <= mPionTpcNSigmaKaon[0]) ||
+           (nSKaon >= mPionTpcNSigmaKaon[1]) ) &&
+          ( (nSProton <= mPionTpcNSigmaProton[0]) ||
+           (nSProton >= mPionTpcNSigmaProton[1]) ) );
 }
 
 //_________________
@@ -511,16 +523,16 @@ Bool_t StFemtoDstInclusiveSelector::IsTpcKaon(StMuTrack *pTrk) {
   Float_t nSPion = pTrk->nSigmaPion();
   Float_t nSKaon = pTrk->nSigmaKaon();
   Float_t nSProton = pTrk->nSigmaProton();
-  
+
   return ( (p >= mKaonTpcMomentum[0]) && (p <= mKaonTpcMomentum[1]) &&
-	   ( (nSKaon >= mKaonTpcNSigmaKaon[0]) &&
-	     (nSKaon <= mKaonTpcNSigmaKaon[1]) ) &&
-	   ( (nSElectron <= mKaonTpcNSigmaElectron[0]) ||
-	     (nSElectron >= mKaonTpcNSigmaElectron[1]) ) &&
-	   ( (nSPion <= mKaonTpcNSigmaPion[0]) ||
-	     (nSPion >= mKaonTpcNSigmaPion[1]) ) &&
-	   ( (nSProton <= mKaonTpcNSigmaProton[0]) ||
-	     (nSProton >= mKaonTpcNSigmaProton[1]) ) );
+          ( (nSKaon >= mKaonTpcNSigmaKaon[0]) &&
+           (nSKaon <= mKaonTpcNSigmaKaon[1]) ) &&
+          ( (nSElectron <= mKaonTpcNSigmaElectron[0]) ||
+           (nSElectron >= mKaonTpcNSigmaElectron[1]) ) &&
+          ( (nSPion <= mKaonTpcNSigmaPion[0]) ||
+           (nSPion >= mKaonTpcNSigmaPion[1]) ) &&
+          ( (nSProton <= mKaonTpcNSigmaProton[0]) ||
+           (nSProton >= mKaonTpcNSigmaProton[1]) ) );
 }
 
 //_________________
@@ -530,37 +542,37 @@ Bool_t StFemtoDstInclusiveSelector::IsTpcProton(StMuTrack *pTrk) {
   Float_t nSPion = pTrk->nSigmaPion();
   Float_t nSKaon = pTrk->nSigmaKaon();
   Float_t nSProton = pTrk->nSigmaProton();
-  
+
   return ( (p >= mProtonTpcMomentum[0]) && (p <= mProtonTpcMomentum[1]) &&
-	   ( (nSProton >= mProtonTpcNSigmaProton[0]) &&
-	     (nSProton <= mProtonTpcNSigmaProton[1]) ) &&
-	   ( (nSElectron <= mProtonTpcNSigmaElectron[0]) ||
-	     (nSElectron >= mProtonTpcNSigmaElectron[1]) ) &&
-	   ( (nSPion <= mProtonTpcNSigmaPion[0]) ||
-	     (nSPion >= mProtonTpcNSigmaPion[1]) ) &&
-	   ( (nSKaon <= mProtonTpcNSigmaKaon[0]) ||
-	     (nSKaon >= mProtonTpcNSigmaKaon[1]) ) );
+          ( (nSProton >= mProtonTpcNSigmaProton[0]) &&
+           (nSProton <= mProtonTpcNSigmaProton[1]) ) &&
+          ( (nSElectron <= mProtonTpcNSigmaElectron[0]) ||
+           (nSElectron >= mProtonTpcNSigmaElectron[1]) ) &&
+          ( (nSPion <= mProtonTpcNSigmaPion[0]) ||
+           (nSPion >= mProtonTpcNSigmaPion[1]) ) &&
+          ( (nSKaon <= mProtonTpcNSigmaKaon[0]) ||
+           (nSKaon >= mProtonTpcNSigmaKaon[1]) ) );
 }
 
 //_________________
 Bool_t StFemtoDstInclusiveSelector::IsTofPion(Float_t mSqr, Float_t p) {
 
   return ( (p >= mPionTofMomentum[0]) && (p <= mPionTofMomentum[1]) &&
-	   (mSqr >= mPionMassSqr[0]) && (mSqr <= mPionMassSqr[1]) );
+          (mSqr >= mPionMassSqr[0]) && (mSqr <= mPionMassSqr[1]) );
 }
 
 //_________________
 Bool_t StFemtoDstInclusiveSelector::IsTofKaon(Float_t mSqr, Float_t p) {
 
   return ( (p >= mKaonTofMomentum[0]) && (p <= mKaonTofMomentum[1]) &&
-	   (mSqr >= mKaonMassSqr[0]) && (mSqr <= mKaonMassSqr[1]) );
+          (mSqr >= mKaonMassSqr[0]) && (mSqr <= mKaonMassSqr[1]) );
 }
 
 //_________________
 Bool_t StFemtoDstInclusiveSelector::IsTofProton(Float_t mSqr, Float_t p) {
 
   return ( (p >= mProtonTofMomentum[0]) && (p <= mProtonTofMomentum[1]) &&
-	   (mSqr >= mProtonMassSqr[0]) && (mSqr <= mProtonMassSqr[1]) );
+          (mSqr >= mProtonMassSqr[0]) && (mSqr <= mProtonMassSqr[1]) );
 }
 
 //_________________

+ 1 - 0
StRoot/StFemtoDstMaker/StFemtoDstInclusiveSelector.h

@@ -207,6 +207,7 @@ class StFemtoDstInclusiveSelector : public StMaker {
   StFemtoTrack *mFTrack;
   
   //Internal class methods
+  Bool_t AcceptRefMult2(StMuTrack *trk);
   Bool_t AcceptTrigger(StMuEvent *muEvent);
   Bool_t AcceptPrimVtx(StThreeVectorF vtxPos, Float_t vpdVz);
   Bool_t AcceptRunId(Int_t runId);   //For StRefMultCorr (isBadRun)

+ 172 - 162
StRoot/StFemtoDstMaker/StFemtoDstMaker.cxx

@@ -6,77 +6,77 @@
 
 ClassImp(StFemtoDstMaker)
 
-// Set maximum file size to 1.9 GB (Root has a 2GB limit)
+  // Set maximum file size to 1.9 GB (Root has a 2GB limit)
 #define MAXFILESIZE 1900000000
 
-//_________________
-StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
-				 const Char_t *oFileName) {
-
-  std::cout << "StFemtoDstMaker::StFemtoDstMaker - Creating an instance...";
-  mOutFileName = oFileName;
-  mMuDstMaker = muMaker;
-  mRefMultCorrUtil = NULL;
-  mCompression = 9;
+  //_________________
+  StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
+                                   const Char_t *oFileName) {
 
-  mCollisionType = false;     //0-pp, 1-AuAu
+    std::cout << "StFemtoDstMaker::StFemtoDstMaker - Creating an instance...";
+    mOutFileName = oFileName;
+    mMuDstMaker = muMaker;
+    mRefMultCorrUtil = NULL;
+    mCompression = 9;
 
-  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] = -10.;
-  mPrimVtxVpdVzDiff[1]= 10.;
-  mPrimVtxXShift = 0.;
-  mPrimVtxYShift = 0.;
-
-  //Initialize single-particle cut variables
-  mTrackMomentum[0] = 0.1;
-  mTrackMomentum[1] = 2.;
-  mTrackDca[0] = 0.;
-  mTrackDca[1] = 5.;
-  mTrackDcaGlobal[0] = 0.;
-  mTrackDcaGlobal[1] = 5.;
-  mTrackNHits[0] = 15;
-  mTrackNHits[1] = 50;
-  mTrackNHitsFit[0] = 15;
-  mTrackNHitsFit[1] = 50;
-
-  mTrackEta[0] = -1.1;
-  mTrackEta[1] = 1.1;
-  mTrackFlag[0] = 0;
-  mTrackFlag[1] = 1000;
-
-  //Flags and cuts for exclusion
-  mRemovePions = false;
-  mRemoveKaons = false;
-  mRemoveProtons = false;
-  mPthresh   = 0.4;
-  mTpcPionNSigma[0] = -1.5;
-  mTpcPionNSigma[1] = 1.5;
-  mTofPionMassSqr[0] = -0.1;
-  mTofPionMassSqr[1] = 0.15;
-  mTpcKaonNSigma[0] = -1.5;
-  mTpcKaonNSigma[1] = 1.5;
-  mTofKaonMassSqr[0] = 0.2;
-  mTofKaonMassSqr[1] = 0.35;
-  mTpcProtonNSigma[0] = -1.5;
-  mTpcProtonNSigma[1] = 1.5;
-  mTofProtonMassSqr[0] = 0.7;
-  mTofProtonMassSqr[1] = 1.1;  
+    mCollisionType = false;     //0-pp, 1-AuAu
 
-  std::cout << "\t[DONE]" << std::endl;
-}
+    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] = -10.;
+    mPrimVtxVpdVzDiff[1]= 10.;
+    mPrimVtxXShift = 0.;
+    mPrimVtxYShift = 0.;
+
+    //Initialize single-particle cut variables
+    mTrackMomentum[0] = 0.1;
+    mTrackMomentum[1] = 2.;
+    mTrackDca[0] = 0.;
+    mTrackDca[1] = 5.;
+    mTrackDcaGlobal[0] = 0.;
+    mTrackDcaGlobal[1] = 5.;
+    mTrackNHits[0] = 15;
+    mTrackNHits[1] = 50;
+    mTrackNHitsFit[0] = 15;
+    mTrackNHitsFit[1] = 50;
+
+    mTrackEta[0] = -1.1;
+    mTrackEta[1] = 1.1;
+    mTrackFlag[0] = 0;
+    mTrackFlag[1] = 1000;
+
+    //Flags and cuts for exclusion
+    mRemovePions = false;
+    mRemoveKaons = false;
+    mRemoveProtons = false;
+    mPthresh   = 0.4;
+    mTpcPionNSigma[0] = -1.5;
+    mTpcPionNSigma[1] = 1.5;
+    mTofPionMassSqr[0] = -0.1;
+    mTofPionMassSqr[1] = 0.15;
+    mTpcKaonNSigma[0] = -1.5;
+    mTpcKaonNSigma[1] = 1.5;
+    mTofKaonMassSqr[0] = 0.2;
+    mTofKaonMassSqr[1] = 0.35;
+    mTpcProtonNSigma[0] = -1.5;
+    mTpcProtonNSigma[1] = 1.5;
+    mTofProtonMassSqr[0] = 0.7;
+    mTofProtonMassSqr[1] = 1.1;  
+
+    std::cout << "\t[DONE]" << std::endl;
+  }
 
 //_________________
 StFemtoDstMaker::~StFemtoDstMaker() {
@@ -87,8 +87,8 @@ StFemtoDstMaker::~StFemtoDstMaker() {
 Int_t StFemtoDstMaker::Init() {
 
   std::cout << "StFemtoDstMaker::Init - Initializing the maker"
-	    << std::endl 
-	    << "Creating output file: " << mOutFileName;
+      << std::endl 
+      << "Creating output file: " << mOutFileName;
 
   mFemtoEvent = new StFemtoEvent();
   mOutFile = new TFile(mOutFileName, "RECREATE");
@@ -110,7 +110,7 @@ Int_t StFemtoDstMaker::Init() {
 
   mNBytes = 0;
   std::cout << "StFemtoDstMaker::Init - Initialization has been finished"
-	    << std::endl;
+      << std::endl;
   return StMaker::Init();
 }
 
@@ -131,7 +131,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;
   }
 
@@ -139,7 +139,7 @@ Int_t StFemtoDstMaker::Make() {
   mMuDst = (StMuDst*)GetInputDS("MuDst");
   if(!mMuDst) {
     gMessMgr->Warning() << "StFemtoDstMaker::Make [WARNING] - No MuDst has been found" 
-			<< endm;
+        << endm;
     return kStOk;
   }
 
@@ -148,7 +148,7 @@ Int_t StFemtoDstMaker::Make() {
   if(!AcceptTrigger(mMuEvent)) { //If trigger is not found
     return kStOk;
   }
-  
+
   //Multiplicity cannot be negative
   if(mMuEvent->refMult() < 0) {
     return kStOk;
@@ -217,13 +217,13 @@ Int_t StFemtoDstMaker::Make() {
     if(mCollisionType == true) { //AuAu collisions
       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());
       }
       mCent9 = mRefMultCorrUtil->getCentralityBin9();
       mCent16 = mRefMultCorrUtil->getCentralityBin16();
@@ -245,43 +245,44 @@ Int_t StFemtoDstMaker::Make() {
     mFemtoEvent->SetRefMultNeg(mMuEvent->refMultNeg());
     //mFemtoEvent->SetZDCe(mMuEvent->zdcTriggerDetector().adc(4));
     //mFemtoEvent->SetZDCw(mMuEvent->zdcTriggerDetector().adc(0));
-    mFemtoEvent->SetBTofTrayMultiplicity( mMuEvent->btofTrayMultiplicity() );
+    mFemtoEvent->SetNumberOfBTofHit( mMuDst->numberOfBTofHit() );
     mFemtoEvent->SetNumberOfPrimaryTracks(mNPrimTracks);
     mFemtoEvent->SetNumberOfGlobalTracks(mNGlobTracks);
     mFemtoEvent->SetMagneticField(mMuEvent->magneticField());
     mFemtoEvent->SetPrimaryVertexPosition(mVertPosition.x(),
-					  mVertPosition.y(),
-					  mVertPosition.z());
+                                          mVertPosition.y(),
+                                          mVertPosition.z());
     mFemtoEvent->SetVpdVz(mVpdVz);
     mFemtoEvent->SetTriggerId(mCurrentTrigger);
     mFemtoEvent->SetPrimaryVertexRanking(mRanking);
 
     StFemtoTrack *mFTrack = NULL;
-  
+    unsigned int refMult2 = 0;
     //Loop over primary tracks
     for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++) {
 
       mPrimTrack = mMuDst->primaryTracks(iTrk);
       mGlobTrack = mMuDst->globalTracks(mPrimTrack->index2Global());
 
+      if( AcceptRefMult2(mPrimTrack) ) refMult2++;
       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;
-	}
+        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.;
+        mBeta = -999.;
+        mMassSqr = -999.;
       }
 
       Bool_t mIsTpcPion = false;
@@ -292,30 +293,30 @@ Int_t StFemtoDstMaker::Make() {
       Bool_t mIsTofProton = false;
 
       if(mRemovePions) {
-	mIsTpcPion = IsTpcPion(mPrimTrack);
-	mIsTofPion = IsTofPion(mGlobTrack);
+        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);
+        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);
+        mIsTpcProton = IsTpcProton(mPrimTrack);
+        mIsTofProton = IsTofProton(mGlobTrack);
 
-	if(mIsTpcProton==true || mIsTofProton==true) {
-	  continue;
-	}
+        if(mIsTpcProton==true || mIsTofProton==true) {
+          continue;
+        }
       }
 
       mFTrack = mFemtoEvent->AddFemtoTrack();
@@ -334,20 +335,20 @@ Int_t StFemtoDstMaker::Make() {
       mFTrack->SetMap0( mGlobTrack->topologyMap().data(0) );
       mFTrack->SetMap1( mGlobTrack->topologyMap().data(1) );
       mFTrack->SetP( mPrimTrack->momentum().x(), 
-		     mPrimTrack->momentum().y(), 
-		     mPrimTrack->momentum().z() );
+                    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() );
+                          mGlobTrack->p().y(),
+                          mGlobTrack->p().z() );
 
       if(mIsTofTrack) {
-	mFTrack->SetBeta(mBeta);
+        mFTrack->SetBeta(mBeta);
       }
       else {
-	mFTrack->SetBeta(-999.);
+        mFTrack->SetBeta(-999.);
       }
 
     } //for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++)
@@ -362,6 +363,15 @@ Int_t StFemtoDstMaker::Make() {
 }
 
 //_________________
+Bool_t StFemtoDstMaker::AcceptRefMult2(StMuTrack *trk) {
+  return (trk->flag() > 0 &&
+          trk->charge() != 0 &&
+          trk->nHitsFit() >= 10 &&
+          fabs(trk->eta()) < 1.0 &&
+          trk->dca().mag() < 3.0);
+}
+
+//_________________
 Bool_t StFemtoDstMaker::AcceptTrigger(StMuEvent *muEvent) {
 
   Bool_t mIsGoodTrigger = false;
@@ -371,9 +381,9 @@ Bool_t StFemtoDstMaker::AcceptTrigger(StMuEvent *muEvent) {
   else {
     for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++) {
       if(muEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) {
-	mIsGoodTrigger = true;
-	mCurrentTrigger = mTriggerIdCollection.at(iTrg);
-	break;
+        mIsGoodTrigger = true;
+        mCurrentTrigger = mTriggerIdCollection.at(iTrg);
+        break;
       }
     } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
   }
@@ -383,21 +393,21 @@ Bool_t StFemtoDstMaker::AcceptTrigger(StMuEvent *muEvent) {
 
 //_________________
 Bool_t StFemtoDstMaker::AcceptPrimVtx(StThreeVectorF vtxPos, 
-				      Float_t vpdVz) {
+                                      Float_t vpdVz) {
   Bool_t mIsGoodVtx = false;
   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;
   if( vtxPos.z() >= mPrimVtxZ[0] &&
-      vtxPos.z() <= mPrimVtxZ[1] &&
-      vtxR >= mPrimVtxR[0] &&
-      vtxR <= mPrimVtxR[1] &&
-      vpdDiff >= mPrimVtxVpdVzDiff[0] &&
-      vpdDiff <= mPrimVtxVpdVzDiff[1]) {
+     vtxPos.z() <= mPrimVtxZ[1] &&
+     vtxR >= mPrimVtxR[0] &&
+     vtxR <= mPrimVtxR[1] &&
+     vpdDiff >= mPrimVtxVpdVzDiff[0] &&
+     vpdDiff <= mPrimVtxVpdVzDiff[1]) {
     mIsGoodVtx = true;
   }
-    
+
   return mIsGoodVtx;
 }
 
@@ -406,21 +416,21 @@ Bool_t StFemtoDstMaker::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) {
 
   Bool_t mIsGoodTrack = false;
   if( trk->vertexIndex() == vtxInd &&
-      trk->momentum().mag() >= mTrackMomentum[0] &&
-      trk->momentum().mag() <= mTrackMomentum[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] &&
-      ((Float_t)trk->nHitsFit()/(Float_t)trk->nHitsPoss()) >= 0.51 &&
-      trk->dca(vtxInd).perp() >= mTrackDca[0] &&
-      trk->dca(vtxInd).perp() <= mTrackDca[1] &&
-      trk->dcaGlobal(vtxInd).perp() >= mTrackDcaGlobal[0] &&
-      trk->dcaGlobal(vtxInd).perp() <= mTrackDcaGlobal[1]) {
+     trk->momentum().mag() >= mTrackMomentum[0] &&
+     trk->momentum().mag() <= mTrackMomentum[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] &&
+     ((Float_t)trk->nHitsFit()/(Float_t)trk->nHitsPoss()) >= 0.51 &&
+     trk->dca(vtxInd).perp() >= mTrackDca[0] &&
+     trk->dca(vtxInd).perp() <= mTrackDca[1] &&
+     trk->dcaGlobal(vtxInd).perp() >= mTrackDcaGlobal[0] &&
+     trk->dcaGlobal(vtxInd).perp() <= mTrackDcaGlobal[1]) {
 
     mIsGoodTrack = true;
   }
@@ -429,18 +439,18 @@ Bool_t StFemtoDstMaker::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) {
   //but if TOF is available, then it should become a mess.
   if(trk->momentum().mag() <= 0.51) {
     if( (trk->nSigmaPion() < -2. ||
-	 trk->nSigmaPion() > 2.) &&
-	TMath::Abs(trk->nSigmaKaon()) > 3. &&
-	TMath::Abs(trk->nSigmaProton()) > 2.)
+         trk->nSigmaPion() > 2.) &&
+       TMath::Abs(trk->nSigmaKaon()) > 3. &&
+       TMath::Abs(trk->nSigmaProton()) > 2.)
       mIsGoodTrack = false;
   }
-      
+
   return mIsGoodTrack;
 }
 
 //_________________
 Bool_t StFemtoDstMaker::IsTofTrack(StMuTrack *trk) { //Only for globals
-  
+
   Bool_t mIsTofHit = false;
   if(trk->btofPidTraits().beta() > 0 &&
      trk->btofPidTraits().timeOfFlight() > 0) {
@@ -458,10 +468,10 @@ Int_t StFemtoDstMaker::Finish() {
   mOutFile->Close();
 
   std::cout << "*************************************" << std::endl
-	    << "StFemtoDstMaker has been finished" << std::endl
-	    << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
-	    << "\t nEventsProcessed: " << mNEventsIn << std::endl
-	    << "*************************************" << std::endl;
+      << "StFemtoDstMaker has been finished" << std::endl
+      << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
+      << "\t nEventsProcessed: " << mNEventsIn << std::endl
+      << "*************************************" << std::endl;
 
   return kStOk;
 }
@@ -478,8 +488,8 @@ Bool_t StFemtoDstMaker::IsTpcPion(StMuTrack *pTrk) {
   Bool_t isPion = false;
 
   if( pTrk->momentum().mag() <= mPthresh && 
-      pTrk->nSigmaPion() >= mTpcPionNSigma[0] &&
-      pTrk->nSigmaPion() <= mTpcPionNSigma[1] ) {
+     pTrk->nSigmaPion() >= mTpcPionNSigma[0] &&
+     pTrk->nSigmaPion() <= mTpcPionNSigma[1] ) {
     isPion = true;
   }
 
@@ -492,8 +502,8 @@ Bool_t StFemtoDstMaker::IsTpcKaon(StMuTrack *pTrk) {
   Bool_t isKaon = false;
 
   if( pTrk->momentum().mag() <= mPthresh && 
-      pTrk->nSigmaKaon() >= mTpcKaonNSigma[0] &&
-      pTrk->nSigmaKaon() <= mTpcKaonNSigma[1] ) {
+     pTrk->nSigmaKaon() >= mTpcKaonNSigma[0] &&
+     pTrk->nSigmaKaon() <= mTpcKaonNSigma[1] ) {
     isKaon = true;
   }
 
@@ -506,8 +516,8 @@ Bool_t StFemtoDstMaker::IsTpcProton(StMuTrack *pTrk) {
   Bool_t isProton = false;
 
   if( pTrk->momentum().mag() <= mPthresh && 
-      pTrk->nSigmaProton() >= mTpcProtonNSigma[0] &&
-      pTrk->nSigmaProton() <= mTpcProtonNSigma[1] ) {
+     pTrk->nSigmaProton() >= mTpcProtonNSigma[0] &&
+     pTrk->nSigmaProton() <= mTpcProtonNSigma[1] ) {
     isProton = true;
   }
 
@@ -521,7 +531,7 @@ Bool_t StFemtoDstMaker::IsTofPion(StMuTrack *gTrk) {
   Float_t beta = -999.;
   Float_t mSqr = -999.;
   if( gTrk->btofPidTraits().beta() > 0 &&
-      gTrk->momentum().mag() > mPthresh ) {
+     gTrk->momentum().mag() > mPthresh ) {
     beta = gTrk->btofPidTraits().beta();
     mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
 
@@ -541,7 +551,7 @@ Bool_t StFemtoDstMaker::IsTofKaon(StMuTrack *gTrk) {
   Float_t beta = -999.;
   Float_t mSqr = -999.;
   if( gTrk->btofPidTraits().beta() > 0 &&
-      gTrk->momentum().mag() > mPthresh ) {
+     gTrk->momentum().mag() > mPthresh ) {
     beta = gTrk->btofPidTraits().beta();
     mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);
 
@@ -561,7 +571,7 @@ Bool_t StFemtoDstMaker::IsTofProton(StMuTrack *gTrk) {
   Float_t beta = -999.;
   Float_t mSqr = -999.;
   if( gTrk->btofPidTraits().beta() > 0 &&
-      gTrk->momentum().mag() > 0.5 ) {
+     gTrk->momentum().mag() > 0.5 ) {
 
     //gTrk->momentum().mag() > mPthresh ) {
     beta = gTrk->btofPidTraits().beta();
@@ -574,10 +584,10 @@ Bool_t StFemtoDstMaker::IsTofProton(StMuTrack *gTrk) {
   }
 
   return isProton;
-}
+  }
 
-//_________________
-void StFemtoDstMaker::SetTriggerId(const UInt_t &id) {
-  mTriggerIdCollection.push_back(id);
-}
+  //_________________
+  void StFemtoDstMaker::SetTriggerId(const UInt_t &id) {
+    mTriggerIdCollection.push_back(id);
+  }
 

+ 1 - 0
StRoot/StFemtoDstMaker/StFemtoDstMaker.h

@@ -161,6 +161,7 @@ public:
   StFemtoEvent *mFemtoEvent;
 
   //Internal class methods
+  Bool_t AcceptRefMult2(StMuTrack *trk);
   Bool_t AcceptTrigger(StMuEvent *muEvent);
   Bool_t AcceptPrimVtx(StThreeVectorF vtxPos, Float_t vpdVz);
   Bool_t AcceptRunId(Int_t runId);   //For StRefMultCorr (isBadRun)

+ 3 - 4
StRoot/StFemtoDstMaker/StFemtoDstPythia6Maker.cxx

@@ -263,11 +263,11 @@ Int_t StFemtoDstPythia6Maker::Make() {
       //Default set to the first bit: primary-vertex-used
       tmap1 |= 1 << 0;
     }
-    else if (ptot < 0.3) { //Tracks may come to the end of the outer sector: R = 2m
+    else if(ptot < 0.3) { //Tracks may come to the end of the outer sector: R = 2m
       tmap1 = 0;
       tmap2 = 0;
-      for (int iBit=0; iBit<32; iBit++)
-      {
+      for(int iBit=0; iBit<32; iBit++) {
+
 	if(randy->Rndm() < 0.95) {
 	  tmap1 |= 1 << iBit;
 	}
@@ -334,7 +334,6 @@ Int_t StFemtoDstPythia6Maker::Make() {
   mFemtoEvent->SetRefMultNeg( pyNegTrkNum );
   //mFemtoEvent->SetZDCe( 45000 );
   //mFemtoEvent->SetZDCw( 45000 );
-  mFemtoEvent->SetBTofTrayMultiplicity( pyRefMult );
   mFemtoEvent->SetNumberOfPrimaryTracks( pyPrimTracksNum );
   mFemtoEvent->SetNumberOfGlobalTracks( pyGlobTracksNum );
   mFemtoEvent->SetMagneticField( mMagneticField );

+ 207 - 167
StRoot/StFemtoDstMaker/StFemtoDstQAMaker.cxx

@@ -5,25 +5,27 @@
 
 ClassImp(StFemtoDstQAMaker)
 
-//_________________
-StFemtoDstQAMaker::StFemtoDstQAMaker(const char *aDir, 
-				     const char *aFileName,
-				     const char *aFilter,
-				     int aMaxFiles) {
-  dir      = string(aDir);
-  fileName = string(aFileName);
-  filter   = aFilter;
-  maxFiles = aMaxFiles;
-  outFileName = "oFemtoQA.root";
-
-  tree       = 0;
-  chain      = 0;
-  femtoEvent = 0;
-  eventId    = 0;
-  mRunIdBins = 7164020;    //Total amount of runs in between Run16 and Run9
-  mRunIdVal[0] = 10016001; //Beginning of Run9
-  mRunIdVal[1] = 17180021; //End of Run16
-}
+  //_________________
+  StFemtoDstQAMaker::StFemtoDstQAMaker(const char *aDir, 
+                                       const char *aFileName,
+                                       const char *aFilter,
+                                       int aMaxFiles) {
+    dir      = string(aDir);
+    fileName = string(aFileName);
+    filter   = aFilter;
+    maxFiles = aMaxFiles;
+    outFileName = "oFemtoQA.root";
+
+    tree       = 0;
+    chain      = 0;
+    femtoEvent = 0;
+    eventId    = 0;
+    mTrig      = 0;
+
+    mRunIdBins = 7164020;    //Total amount of runs in between Run16 and Run9
+    mRunIdVal[0] = 10016001; //Beginning of Run9
+    mRunIdVal[1] = 17180021; //End of Run16
+  }
 
 //_________________
 Int_t StFemtoDstQAMaker::Init() {
@@ -41,15 +43,15 @@ Int_t StFemtoDstQAMaker::Init() {
   Int_t mPrimVtxZBins = 400;
   Double_t mPrimVtxZLo = -200.;
   Double_t mPrimVtxZHi = 200.;
-  Int_t mPrimVertRBins = 80;
+  Int_t mPrimVertRBins = 60;
   Double_t mPrimVertRLo = -0.5;
-  Double_t mPrimVertRHi = 1.5;
+  Double_t mPrimVertRHi = 5.5;
   Int_t mPrimVertXBins = 100;
-  Double_t mPrimVertXLo = -1.;
-  Double_t mPrimVertXHi = 1.;
+  Double_t mPrimVertXLo = -5.;
+  Double_t mPrimVertXHi = 5.;
   Int_t mPrimVertYBins = 100;
-  Double_t mPrimVertYLo = -1.;
-  Double_t mPrimVertYHi = 1.;
+  Double_t mPrimVertYLo = -5.;
+  Double_t mPrimVertYHi = 5.;
   Int_t mNSigmaBins = 200;
   Double_t mNSigmaLo = -10.;
   Double_t mNSigmaHi = 10.;
@@ -82,103 +84,114 @@ Int_t StFemtoDstQAMaker::Init() {
   Double_t mBTofMultHi = 699.5;
 
   // Event histograms
-  hRefMultWeight = new TH1F("hRefMultWeight","RefMultCorrWeight;weight", 
-			    mRefMultWeightBins, mRefMultWeightLo, mRefMultWeightHi);
+  hRefMultWeight = new TH1F("hRefMultWeight",";RefMultCorrWeight;#",
+                            mRefMultWeightBins, mRefMultWeightLo, mRefMultWeightHi);
   hRefMult = new TH1F("hRefMult",
-		      "refMult;RefMult;dN/dRefMult",
-		      mRefMultBins, mRefMultLo, mRefMultHi);
-  hBTofTrayMultVsRefMult = new TH2F("hBTofTrayMultVsRefMult","BTofTrayMult vs. refMult;refMult;BTofTrayMult",
-				    mRefMultBins, mRefMultLo, mRefMultHi,
-				    mBTofMultBins, mBTofMultLo, mBTofMultHi);
-  hRefMultCorr = new TH1F("hRefMultCorr", "RefMultCorr;RefMultCorr",
-			  mRefMultBins, mRefMultLo, mRefMultHi);
+                      ";RefMult;#",
+                      mRefMultBins, mRefMultLo, mRefMultHi);
+  hRefMult2 = new TH1F("hRefMult2",
+                       "|#eta| < 1;RefMult2;#",
+                       mRefMultBins, mRefMultLo, mRefMultHi);
+  hRefMultCorr = new TH1F("hRefMultCorr", ";RefMultCorr;#",
+                          mRefMultBins, mRefMultLo, mRefMultHi);
+  hTofMult = new TH1F("hTofMult", ";TofMult;#",
+                      mRefMultBins, mRefMultLo, mRefMultHi);
+  hTofVsRef = new TH2F("hTofVsRef", ";TofMult;RefMult",
+                       mRefMultBins, mRefMultLo, mRefMultHi,
+                       mRefMultBins, mRefMultLo, mRefMultHi);
+  hTofVsRef2 = new TH2F("hTofVsRef2", ";TofMult;RefMult2",
+                        mRefMultBins, mRefMultLo, mRefMultHi,
+                        mRefMultBins, mRefMultLo, mRefMultHi);
+  hRefVsRef2 = new TH2F("hRefVsRef2", ";RefMult;RefMult2",
+                        mRefMultBins, mRefMultLo, mRefMultHi,
+                        mRefMultBins, mRefMultLo, mRefMultHi);
   hPrimVertZ = new TH1F("hPrimVertZ",
-			"Z coordinate of primary vertex;z (cm);dN/dz",
-			mPrimVtxZBins, mPrimVtxZLo, mPrimVtxZHi);
-  hPrimVertR = new TH1F("hPrimVertR","#sqrt{V_{x}^{2}+V_{y}^{2}} of primary vertex;R (cm)",
-		        mPrimVertRBins, mPrimVertRLo, mPrimVertRHi);
+                        "Z coordinate of primary vertex;z (cm);#",
+                        mPrimVtxZBins, mPrimVtxZLo, mPrimVtxZHi);
+  hPrimVertR = new TH1F("hPrimVertR","#sqrt{V_{x}^{2}+V_{y}^{2}} of primary vertex;R (cm);#",
+                        mPrimVertRBins, mPrimVertRLo, mPrimVertRHi);
   hPrimVertXvsY = new TH2F("hPrimVertXvsY","X vs Y for primary vertex;x (cm);y (cm)",
-			   mPrimVertXBins, mPrimVertXLo, mPrimVertXHi,
-			   mPrimVertYBins, mPrimVertYLo, mPrimVertYHi);
-  hPrimVertVpdVzDiff = new TH1F("hPrimVertVpdVzDiff","z_{TPC} - z_{VPD};z_{TPC} - z_{VPD} (cm)",
-				300, -15., 15.);
+                           mPrimVertXBins, mPrimVertXLo, mPrimVertXHi,
+                           mPrimVertYBins, mPrimVertYLo, mPrimVertYHi);
+  hPrimVertVpdVzDiff = new TH1F("hPrimVertVpdVzDiff","z_{TPC} - z_{VPD};z_{TPC} - z_{VPD} (cm);#",
+                                300, -15., 15.);
   hRefMultVsRunNumer = new TProfile("hRefMultVsRunNumer","refMult vs. Run Number;Run Number;<refMult>",
-				    mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
-  hBTofTrayMultVsRunNumber = new TProfile("hBTofTrayMultVsRunNumber","BTofTrayMult vs. Run Number;Run Number;<btofTrayMult>",
-					  mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
-  hPrimTrackNumVsRunNumber = new TProfile("hPrimTrackNumVsRunNumber","# pTrk vs. Run Number;Run Number;<PrimTrackNum>",
-					  mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
-  hGlobTrackNumVsRunNumber = new TProfile("hGlobTrackNumVsRunNumber","#gTrk vs. Run Number;Run Number;<GlobTrackNum>",
-					  mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
+                                    mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
+  hPrimTrackNumVsRunNumber = new TProfile("hPrimTrackNumVsRunNumber","pTrk vs. Run Number;Run Number;<PrimTrackNum>",
+                                          mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
+  hGlobTrackNumVsRunNumber = new TProfile("hGlobTrackNumVsRunNumber","gTrk vs. Run Number;Run Number;<GlobTrackNum>",
+                                          mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
+  hTofMultVsRunNumber = new TProfile("hTofMultVsRunNumber","TofMult vs. Run Number;Run Number;<TofMult>",
+                                     mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
   // Track histograms
   hNSigmaPionVsMom = new TH2F("hNSigmaPionVsMom", "N#sigma(#pi) vs momenta;p (GeV/c);N#sigma(pi)",
-			      mMomBins,mMomLo,mMomHi,
-			      mNSigmaBins, mNSigmaLo, mNSigmaHi);
+                              mMomBins,mMomLo,mMomHi,
+                              mNSigmaBins, mNSigmaLo, mNSigmaHi);
   hNSigmaKaonVsMom = new TH2F("hNSigmaKaonVsMom", "N#sigma(K) vs momenta;p (GeV/c);N#sigma(K)",
-			      mMomBins,mMomLo,mMomHi,
-			      mNSigmaBins, mNSigmaLo, mNSigmaHi);
+                              mMomBins,mMomLo,mMomHi,
+                              mNSigmaBins, mNSigmaLo, mNSigmaHi);
   hNSigmaProtonVsMom = new TH2F("hNSigmaProtonVsMom", "N#sigma(p) vs momenta;p (GeV/c);N#sigma(p)",
-			      mMomBins,mMomLo,mMomHi,
-			      mNSigmaBins, mNSigmaLo, mNSigmaHi);
+                                mMomBins,mMomLo,mMomHi,
+                                mNSigmaBins, mNSigmaLo, mNSigmaHi);
   hNSigmaElectronVsMom = new TH2F("hNSigmaElectronVsMom", "N#sigma(e) vs momenta;p (GeV/c);N#sigma(e)",
-			      mMomBins,mMomLo,mMomHi,
-			      mNSigmaBins, mNSigmaLo, mNSigmaHi);
+                                  mMomBins,mMomLo,mMomHi,
+                                  mNSigmaBins, mNSigmaLo, mNSigmaHi);
 
-  hP = new TH1F("hP","Momentum;p (GeV/c);dN/dp",mMomBins,mMomLo,mMomHi);
-  hPt = new TH1F("hPt","Transverse momentum;p_{T} (GeV/c);dN/dp_{T}",mMomBins,mMomLo,mMomHi);
+  hP = new TH1F("hP","Momentum;p (GeV/c);#",mMomBins,mMomLo,mMomHi);
+  hPt = new TH1F("hPt","Transverse momentum;p_{T} (GeV/c);#",mMomBins,mMomLo,mMomHi);
 
   hGlobPx = new TH1F("hGlobPx",
-		     "global track p_{x};p_{x} (GeV/c);N",
-		     2*mMomBins,-mMomHi,mMomHi);
+                     "global track p_{x};p_{x} (GeV/c);#",
+                     2*mMomBins,-mMomHi,mMomHi);
   hGlobPy = new TH1F("hGlobPy",
-		     "global track p_{y};p_{y} (GeV/c);N",
-		     2*mMomBins,-mMomHi,mMomHi);
+                     "global track p_{y};p_{y} (GeV/c);#",
+                     2*mMomBins,-mMomHi,mMomHi);
   hGlobPz = new TH1F("hGlobPz",
-		     "global track p_{z};p_{z} (GeV/c);N",
-		     2*mMomBins,-mMomHi,mMomHi);
+                     "global track p_{z};p_{z} (GeV/c);#",
+                     2*mMomBins,-mMomHi,mMomHi);
 
   hGlobX = new TH1F("hGlobX",
-		    "global track first point x-coordinate;x (cm);N",
-		    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
+                    "global track first point x-coordinate;x (cm);#",
+                    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
   hGlobY = new TH1F("hGlobY",
-		    "global track first point y-coordinate;y (cm);N",
-		    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
+                    "global track first point y-coordinate;y (cm);#",
+                    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
   hGlobZ = new TH1F("hGlobZ",
-		    "global track first point z-coordinate;z (cm);N",
-		    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
+                    "global track first point z-coordinate;z (cm);#",
+                    mGlobCoordBins, mGlobCoordLo, mGlobCoordHi);
+
 
-  
   hdEdx = new TH1F("hdEdx", "dE/dx;dE/dx (keV/cm);N", mDedxBins, mDedxLo, mDedxHi);
   hdEdx_p = new TH2F("hdEdx_p", "dE/dx vs momentum;q*p (GeV/c); dE/dx (keV/cm);", 
                      2*mMomBins, -mMomHi, mMomHi, 
-		     mDedxBins, mDedxLo, mDedxHi);
+                     mDedxBins, mDedxLo, mDedxHi);
   hDCAxy = new TH1F("hDCAxy",
-		    "Distance of closest approach (XY plane) for primary track;DCA_{xy} (cm);N",
-		    mDCAxyBins, mDCAxyLo, mDCAxyHi);
+                    "Distance of closest approach (XY plane) for primary track;DCA_{xy} (cm);#",
+                    mDCAxyBins, mDCAxyLo, mDCAxyHi);
   hDCAz = new TH1F("hDCAz",
-		   "Distance of closest approach (Z axis) for primary track;DCA_{z} (cm);N",
-		   mDCAzBins, mDCAzLo, mDCAzHi);
+                   "Distance of closest approach (Z axis) for primary track;DCA_{z} (cm);#",
+                   mDCAzBins, mDCAzLo, mDCAzHi);
   hDCAxyGlob = new TH1F("hDCAxyGlob",
-			"Distance of closest approach (XY plane) for global track;DCA_{xy} (cm);N",
-			mDCAxyBins, mDCAxyLo, mDCAxyHi);
+                        "Distance of closest approach (XY plane) for global track;DCA_{xy} (cm);#",
+                        mDCAxyBins, mDCAxyLo, mDCAxyHi);
   hDCAzGlob = new TH1F("hDCAzGlob",
-		       "Distance of closest approach (Z axis) for global track;DCA_{z} (cm);N",
-		       mDCAzBins, mDCAzLo, mDCAzHi);
+                       "Distance of closest approach (Z axis) for global track;DCA_{z} (cm);#",
+                       mDCAzBins, mDCAzLo, mDCAzHi);
   hMSqr_P = new TH2F("hMSqr_P","m^{2} vs p;p (GeV/c);m^{2} (GeV/c^{2})^{2}",
                      2*mMomBins, -mMomHi, mMomHi, 
-		     mMSqrBins, mMSqrLo, mMSqrHi);
+                     mMSqrBins, mMSqrLo, mMSqrHi);
   hInvBeta_P = new TH2F("hInvBeta_P",
-			"1/#beta vs p*q;q*p (GeV/c);1/#beta",
-			2*mMomBins, -mMomHi, mMomHi, 
-			mInvBetaBins, mInvBetaLo, mInvBetaHi);
+                        "1/#beta vs p*q;q*p (GeV/c);1/#beta",
+                        2*mMomBins, -mMomHi, mMomHi, 
+                        mInvBetaBins, mInvBetaLo, mInvBetaHi);
   hBeta_P = new TH2F("hBeta_P",
-		     "#beta vs momentum;q*p (GeV/c);#beta",
-		     2*mMomBins, -mMomHi, mMomHi, 
-		     mBetaBins, mBetaLo, mBetaHi);
+                     "#beta vs momentum;q*p (GeV/c);#beta",
+                     2*mMomBins, -mMomHi, mMomHi, 
+                     mBetaBins, mBetaLo, mBetaHi);
   hEta = new TH1F("hEta",
-		  "#eta;#eta;N",
-		  100, -1., 1.);
-  
+                  "#eta;#eta;N",
+                  100, -1., 1.);
+
   return kStOk;
 }
 
@@ -194,7 +207,7 @@ Int_t StFemtoDstQAMaker::Finish() {
 
 //_________________
 int StFemtoDstQAMaker::FillChain(TChain *aChain, char *aDir,
-				 const char *aFilter, int aMaxFiles) {
+                                 const char *aFilter, int aMaxFiles) {
 
   TSystem *gSystem = new TSystem();
   if (!gSystem) {
@@ -204,11 +217,11 @@ int StFemtoDstQAMaker::FillChain(TChain *aChain, char *aDir,
   void *lDir = gSystem->OpenDirectory(aDir);
   if (!lDir) {
     std::cout << "StFemtoDstQAMaker[ERROR]: can't open directory " 
-	      << aDir << std::endl;
+        << aDir << std::endl;
     delete gSystem;
     return 0;
   }
-  
+
   int lCount = 0;
   char *lFileName;
 
@@ -218,33 +231,33 @@ int StFemtoDstQAMaker::FillChain(TChain *aChain, char *aDir,
     if (strstr(lFileName, aFilter)) {
       char *lFullName = gSystem->ConcatFileName(aDir, lFileName);
       std::cout << "StFemtoDstQAMaker[INFO]: Adding "
-		<< lFullName << " to the chain" << std::endl;
+          << lFullName << " to the chain" << std::endl;
       aChain->Add(lFullName);
       delete lFullName;
       lCount++;
       if ((aMaxFiles > 0) && (lCount > aMaxFiles)) break;
     }
   }
-  
+
   std::cout << "StFemtoDstQAMaker[INFO]: Added " << lCount 
-	    << " files to the chain" << std::endl;
+      << " files to the chain" << std::endl;
   delete gSystem;
   return lCount;
 }
 
 //_________________
 int StFemtoDstQAMaker::FillChain(TChain *aChain, const char *aFileName,
-				 int aMaxFiles) {
+                                 int aMaxFiles) {
   std::ifstream lInStream(aFileName);
 
 
   if (!lInStream.is_open()) {
     std::cout << "StFemtoDstQAMaker[ERROR]: can't open file list: "
-	      << aFileName << std::endl;
+        << aFileName << std::endl;
     return 0;
   }
   std::cout << "StFemtoDstQAMaker[INFO]: using file list: "
-	    << aFileName << std::endl;
+      << aFileName << std::endl;
 
   int lCount = 0;
   char buf[0xFF];
@@ -259,26 +272,26 @@ int StFemtoDstQAMaker::FillChain(TChain *aChain, const char *aFileName,
   }
 
   std::cout << "StFemtoDstQAMaker[INFO]: Added " << lCount
-	    << " files to the chain" << std::endl;
+      << " files to the chain" << std::endl;
   return lCount;
 }
 
 //_________________
 int StFemtoDstQAMaker::InitRead(string aDir, string aFileName,
-				string aFilter, int aMaxFiles) {
+                                string aFilter, int aMaxFiles) {
   eventId = 0;
   //  femtoEvent = new StFemtoEvent();
   chain = new TChain("StFemtoDst", "StFemtoDst");
-  
+
   std::cout << "StFemtoDstQAMaker[INFO]: Call InitRead" << std::endl;
-  
+
   int lNumFiles = 0;
   if (!aFileName.empty()) {
     if (strstr(aFileName.c_str(), ".lis") ||
-	strstr(aFileName.c_str(), ".list")) {
+        strstr(aFileName.c_str(), ".list")) {
       lNumFiles = FillChain(chain, 
-			    (aDir + aFileName).c_str(), 
-			    aMaxFiles);
+                            (aDir + aFileName).c_str(), 
+                            aMaxFiles);
     }
     else {
       chain->Add((aDir + aFileName).c_str());
@@ -286,7 +299,7 @@ int StFemtoDstQAMaker::InitRead(string aDir, string aFileName,
     }
   } else {
     lNumFiles = FillChain(chain, (char *)aDir.c_str(), // not cool
-			  aFilter.c_str(), aMaxFiles);
+                          aFilter.c_str(), aMaxFiles);
   }
   chain->SetBranchAddress("StFemtoEvent", &femtoEvent);
   nEvents = (unsigned int)chain->GetEntries();
@@ -302,7 +315,7 @@ void StFemtoDstQAMaker::UninitRead() {
 
 //_________________
 void StFemtoDstQAMaker::AddTrigId(unsigned int aTrigId) {
-  trigColl.push_back(aTrigId);
+  mTrigIdCol.push_back(aTrigId);
 }
 
 //_________________
@@ -334,84 +347,111 @@ Int_t StFemtoDstQAMaker::Make() {
   }
 
   if (femtoEvent) {
-    if (!trigColl.empty()) {
-      bool lTrigPass = false;
+    if (!mTrigIdCol.empty()) {
       unsigned int lTrigId = femtoEvent->GetTriggerId();
-      
-      for (unsigned int i = 0; i < trigColl.size(); i++) {
-	if (lTrigId == trigColl[i]) {
-	  lTrigPass = true;
-	  break;
-	}
-      }
-      if (!lTrigPass) return 0;
+      if (((lTrigId & mTrig) == 0) || ((lTrigId & (!mTrig)) != 0))
+        return 0;
     }
 
     // fill histograms
     hRefMultWeight->Fill(femtoEvent->GetRefMultCorrWeight());
     hRefMult->Fill(femtoEvent->GetRefMult());
+    hRefMult2->Fill(femtoEvent->GetRefMult2());
     hRefMultCorr->Fill(femtoEvent->GetRefMultCorr());
+    hTofMult->Fill(femtoEvent->GetNumberOfBTofHit());
+    hTofVsRef->Fill(femtoEvent->GetNumberOfBTofHit(),
+                    femtoEvent->GetRefMult());
+    hTofVsRef2->Fill(femtoEvent->GetNumberOfBTofHit(),
+                     femtoEvent->GetRefMult2());
+    hRefVsRef2->Fill(femtoEvent->GetRefMult(),
+                     femtoEvent->GetRefMult2());
     hPrimVertZ->Fill(femtoEvent->GetVertexPositionZ());
     hPrimVertR->Fill(femtoEvent->GetVertexPositionR());
     hPrimVertXvsY->Fill(femtoEvent->GetVertexPositionX(),
-			femtoEvent->GetVertexPositionY());
+                        femtoEvent->GetVertexPositionY());
     hPrimVertVpdVzDiff->Fill(femtoEvent->GetVertexPositionZ() - femtoEvent->GetVpdVz());
     hRefMultVsRunNumer->Fill(femtoEvent->GetRunId(),femtoEvent->GetRefMult());
-    hBTofTrayMultVsRunNumber->Fill(femtoEvent->GetRunId(),femtoEvent->GetBTofTrayMultiplicity());
-    hPrimTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),femtoEvent->GetNumberOfPrimaryTracks());
-    hGlobTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),femtoEvent->GetNumberOfGlobalTracks());
-    hBTofTrayMultVsRefMult->Fill(femtoEvent->GetRefMult(),femtoEvent->GetBTofTrayMultiplicity());
+    hPrimTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),
+                                   femtoEvent->GetNumberOfPrimaryTracks());
+    hGlobTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),
+                                   femtoEvent->GetNumberOfGlobalTracks());
+    hTofMultVsRunNumber->Fill(femtoEvent->GetRunId(),
+                              femtoEvent->GetNumberOfBTofHit());
 
     TClonesArray *lTracks = femtoEvent->GetFemtoTracks();
 
     if (lTracks) {
       int lNTracks = lTracks->GetEntries();
       for (int i = 0; i < lNTracks; i++) {
-	StFemtoTrack *lTrack = (StFemtoTrack *)lTracks->At(i);
-
-	hNSigmaPionVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaPion());
-	hNSigmaKaonVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaKaon());
-	hNSigmaProtonVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaProton());
-	hNSigmaElectronVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaElectron());
-
-	hGlobPx->Fill( lTrack->GetGlobalPx() );
-	hGlobPy->Fill( lTrack->GetGlobalPy() );
-	hGlobPz->Fill( lTrack->GetGlobalPz() );
-	hGlobX->Fill( femtoEvent->GetVertexPositionX() + lTrack->GetDCAx() );
-	hGlobY->Fill( femtoEvent->GetVertexPositionY() + lTrack->GetDCAy() );
-	hGlobZ->Fill( femtoEvent->GetVertexPositionZ() + lTrack->GetDCAz() );
-
-	StThreeVectorF vMomGlob( lTrack->GetGlobalPx(),
-				 lTrack->GetGlobalPy(),
-				 lTrack->GetGlobalPz() );
-	StThreeVectorF vOrigGlob( femtoEvent->GetVertexPositionX() + lTrack->GetDCAx(),
-				  femtoEvent->GetVertexPositionY() + lTrack->GetDCAy(),
-				  femtoEvent->GetVertexPositionZ() + lTrack->GetDCAz() );
-	StPhysicalHelixD helixGlob( vMomGlob, vOrigGlob,
-				    femtoEvent->GetMagneticField()*kilogauss,
-				    lTrack->GetCharge() );
-	
-	hdEdx->Fill( lTrack->GetDedx()*1e6 );
-	hdEdx_p->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetDedx()*1e6 );
-	hDCAxy->Fill( lTrack->GetDCAxy() );
-	hDCAz->Fill( lTrack->GetDCAz() );
-	hDCAxyGlob->Fill( lTrack->GetDCAxyGlobal() );
-	hDCAzGlob->Fill( lTrack->GetDCAzGlobal() );
-	hEta->Fill( lTrack->GetEta() );
-	hP->Fill (lTrack->GetPtot());
-	hPt->Fill (lTrack->GetPt());
-	
-
-	if (lTrack->GetBeta() > 0.) {
-	  Float_t mInvBeta = 1./lTrack->GetBeta();
-	  
-	  hMSqr_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetMassSqr() );
-	  hInvBeta_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), mInvBeta );
-	  hBeta_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetBeta() );
-	}
-		
+        StFemtoTrack *lTrack = (StFemtoTrack *)lTracks->At(i);
+
+        hNSigmaPionVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaPion());
+        hNSigmaKaonVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaKaon());
+        hNSigmaProtonVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaProton());
+        hNSigmaElectronVsMom->Fill(lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetNSigmaElectron());
+
+        hGlobPx->Fill( lTrack->GetGlobalPx() );
+        hGlobPy->Fill( lTrack->GetGlobalPy() );
+        hGlobPz->Fill( lTrack->GetGlobalPz() );
+        hGlobX->Fill( femtoEvent->GetVertexPositionX() + lTrack->GetDCAx() );
+        hGlobY->Fill( femtoEvent->GetVertexPositionY() + lTrack->GetDCAy() );
+        hGlobZ->Fill( femtoEvent->GetVertexPositionZ() + lTrack->GetDCAz() );
+
+        StThreeVectorF vMomGlob( lTrack->GetGlobalPx(),
+                                lTrack->GetGlobalPy(),
+                                lTrack->GetGlobalPz() );
+        StThreeVectorF vOrigGlob( femtoEvent->GetVertexPositionX() + lTrack->GetDCAx(),
+                                 femtoEvent->GetVertexPositionY() + lTrack->GetDCAy(),
+                                 femtoEvent->GetVertexPositionZ() + lTrack->GetDCAz() );
+        StPhysicalHelixD helixGlob( vMomGlob, vOrigGlob,
+                                   femtoEvent->GetMagneticField()*kilogauss,
+                                   lTrack->GetCharge() );
+
+        hdEdx->Fill( lTrack->GetDedx()*1e6 );
+        hdEdx_p->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetDedx()*1e6 );
+        hDCAxy->Fill( lTrack->GetDCAxy() );
+        hDCAz->Fill( lTrack->GetDCAz() );
+        hDCAxyGlob->Fill( lTrack->GetDCAxyGlobal() );
+        hDCAzGlob->Fill( lTrack->GetDCAzGlobal() );
+        hEta->Fill( lTrack->GetEta() );
+
+        if (lTrack->GetBeta() > 0.) {
+          Float_t mInvBeta = 1./lTrack->GetBeta();
+
+          hMSqr_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetMassSqr() );
+          hInvBeta_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), mInvBeta );
+          hBeta_P->Fill( lTrack->GetCharge()*lTrack->GetPtot(), lTrack->GetBeta() );
+        }
+
       }
     }
   }
   return kStOk;
 }
+
+  
+void StFemtoDstQAMaker::SetTrigger(unsigned int trig, char mode)
+{
+  unsigned int trigColSize = mTrigIdCol.size();
+  for (unsigned int iTrig = 0; iTrig < trigColSize; iTrig++) {
+    if (mTrigIdCol[iTrig] == trig) {
+      unsigned int mask = 1 << iTrig;
+      if (!mode)
+        mTrig = mTrig & (!mask);
+      else
+        mTrig = mTrig | mask;
+      break;
+    }
+  }
+
+  //
+  // Output trigger mask
+  //
+  char bitMask[33] = {0};
+  for (unsigned int iBit = 0; iBit < 32; iBit++) {
+    bitMask[31 - iBit] = '0' + ((mTrig & (1 << iBit)) >> iBit);
+  }
+  std::cout << "StFemtoDstQAMaker[INFO]: Trigger mask = "
+            << mTrig << " bit = "
+            << bitMask << std::endl;
+}

+ 19 - 7
StRoot/StFemtoDstMaker/StFemtoDstQAMaker.h

@@ -30,6 +30,7 @@ class StFemtoDstQAMaker : public StMaker {
   void SetOutFileName(const char *oFileName) {outFileName = oFileName;}
   //Should be set for each year and run period. The values can be found in: online.star.bnl.gov->RunLog Browser
   void SetRunIdParameters(unsigned int nbins, double lo, double hi) { mRunIdBins=nbins; mRunIdVal[0]=lo; mRunIdVal[1]=hi; }
+  void SetTrigger(unsigned int trig, char mode);
 
  private:
   string dir;
@@ -46,21 +47,31 @@ class StFemtoDstQAMaker : public StMaker {
   unsigned int mRunIdBins;
   double mRunIdVal[2];
 
-  //Event histograms
+  //
+  // Event histograms
+  //
   TH1F *hRefMultCorr;
   TH1F *hRefMult;
+  TH1F *hRefMult2;
   TH1F *hRefMultWeight;
+  TH1F *hTofMult;
+  TH2F *hTofVsRef;
+  TH2F *hTofVsRef2;
+  TH2F *hRefVsRef2;
   TH1F *hPrimVertZ;
   TH1F *hPrimVertR;
   TH2F *hPrimVertXvsY;
   TH1F *hPrimVertVpdVzDiff;
-  TH2F *hBTofTrayMultVsRefMult;
-  //Here is the way of how to look on the quantities over the run
+  //
+  // Here is the way of how to look on the quantities over the run
+  //
   TProfile *hRefMultVsRunNumer;
-  TProfile *hBTofTrayMultVsRunNumber;
   TProfile *hPrimTrackNumVsRunNumber;
   TProfile *hGlobTrackNumVsRunNumber;
-  //Track histograms
+  TProfile *hTofMultVsRunNumber;
+  //
+  // Track histograms
+  //
   TH2F *hNSigmaPionVsMom;
   TH2F *hNSigmaKaonVsMom;
   TH2F *hNSigmaProtonVsMom;
@@ -88,7 +99,8 @@ class StFemtoDstQAMaker : public StMaker {
   TH2F *hBeta_P;
   TH1F *hEta;
 
-  std::vector<unsigned int> trigColl;
+  std::vector<unsigned int> mTrigIdCol;
+  unsigned int              mTrig;
 
   int maxFiles;
   
@@ -100,7 +112,7 @@ class StFemtoDstQAMaker : public StMaker {
 		const char *aFilter, int aMaxFiles);
   void UninitRead();
   
-  ClassDef(StFemtoDstQAMaker, 2)
+  ClassDef(StFemtoDstQAMaker, 3)
 };
 
 #endif

+ 11 - 0
StRoot/StFemtoDstMaker/StFemtoEvent.cxx

@@ -85,6 +85,17 @@ void StFemtoEvent::SetRefMult(Int_t mult) {
 }
 
 //_________________
+void StFemtoEvent::SetRefMult2(Int_t mult) {
+
+  if(mult<0 || mult>60000) {
+    mRefMult2 = USHORTMAX;
+  }
+  else {
+    mRefMult2 = (UShort_t)mult;
+  }
+}
+
+//_________________
 void StFemtoEvent::SetRefMultCorr(Float_t rf) {
 
   if(rf<0. || rf >6000) {

+ 8 - 4
StRoot/StFemtoDstMaker/StFemtoEvent.h

@@ -27,11 +27,13 @@ class StFemtoEvent : public TObject {
   void SetRefMult(Int_t mult);
   void SetRefMultCorr(Float_t rf);
   void SetRefMultCorrWeight(Float_t w);
+  void SetRefMult2(Int_t mult);
   void SetCent9(Int_t c)               { mCent9 = (Char_t)c; }
   void SetCent16(Int_t c)              { mCent16 = (Char_t)c; }
   void SetRefMultPos(Int_t m)          { mRefMultPos = (UShort_t)m; }
   void SetRefMultNeg(Int_t m)          { mRefMultNeg = (UShort_t)m; }
-  void SetBTofTrayMultiplicity(Int_t m)  { mBTofMultiplicity = (UShort_t)m; }
+  void SetNumberOfBTofHit(Int_t m)     { if(m<=0) {mNumberOfBTofHit=0;} else {mNumberOfBTofHit = (UShort_t)m;} }
+  void SetNumberOfBTofHit(UInt_t m)    { mNumberOfBTofHit = (UShort_t)m; }
   void SetNumberOfPrimaryTracks(Int_t m) { mNumberOfPrimaryTracks = (UShort_t)m; }
   void SetNumberOfGlobalTracks(Int_t m)  { mNumberOfGlobalTracks = (UShort_t)m; }
   void SetMagneticField(Float_t f);
@@ -52,13 +54,14 @@ class StFemtoEvent : public TObject {
   UInt_t   GetRefMult() const          { return (UInt_t)mRefMult; }
   Float_t  GetRefMultCorr() const      { return (Float_t)mRefMultCorr * 0.1; }
   Float_t  GetRefMultCorrWeight() const{ return (Float_t)mRefMultCorrWeight * 0.0001; }
+  UInt_t   GetRefMult2() const         { return (UInt_t)mRefMult2; }
   UShort_t GetCent9() const            { return (UShort_t)mCent9; }
   UShort_t GetCent16() const           { return (UShort_t)mCent16; }
   UShort_t GetRefMultPos() const       { return mRefMultPos; }
   UShort_t GetRefMultNeg() const       { return mRefMultNeg; }
   Float_t  GetZDCe() const             { return (Float_t)30000; }
   Float_t  GetZDCw() const             { return (Float_t)30000; }
-  UShort_t GetBTofTrayMultiplicity() const  { return mBTofMultiplicity; }
+  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; }
@@ -81,13 +84,14 @@ class StFemtoEvent : public TObject {
   UShort_t mRefMult;             //mRefMult
   UShort_t mRefMultCorr;         //mRefMultCorr * 10
   UShort_t mRefMultCorrWeight;   //mRefMultCorrWight * 10000
+  UShort_t mRefMult2;            //mRefMult2 for |eta| < 1.0
   Char_t   mCent9;
   Char_t   mCent16;
   UShort_t mRefMultPos;
   UShort_t mRefMultNeg;
   UShort_t mNumberOfPrimaryTracks;
   UShort_t mNumberOfGlobalTracks;
-  UShort_t mBTofMultiplicity;   
+  UShort_t mNumberOfBTofHit;   
   Short_t  mMagField;               //Field*1000
   Float_t  mVertexPositionX;        //VtxX
   Float_t  mVertexPositionY;        //VtxY
@@ -99,7 +103,7 @@ class StFemtoEvent : public TObject {
   UShort_t mNFemtoTracks;
   TClonesArray *mFemtoTracks;
 
-  ClassDef(StFemtoEvent, 2)
+  ClassDef(StFemtoEvent, 4)
 };
 
 #endif

+ 1 - 1
StRoot/StFemtoDstMaker/StFemtoTrack.cxx

@@ -92,7 +92,7 @@ Float_t StFemtoTrack::GetPtot() const {
 Float_t StFemtoTrack::GetEta() const {
 
   float ptot = TMath::Sqrt(mMomX*mMomX + mMomY*mMomY + mMomZ*mMomZ);
-  float eta = 0.5*TMath::Log( (ptot + mMomZ) / (ptot - mMomZ));
+  float eta = 0.5*TMath::Log( (ptot - mMomZ) / (ptot + mMomZ));
   return eta;
 }
 

+ 3 - 3
StRoot/StFemtoDstMaker/StFemtoTrack.h

@@ -56,7 +56,7 @@ class StFemtoTrack : public TObject {
 
   //Getters
   Short_t   GetId() const                    { return (Short_t)mId; }
-  Short_t   GetType() const                  { return (Short_t)1;}
+  Short_t   GetType() const                  { return (Short_t)((mMomX == 0.0 && mMomY == 0.0 && mMomZ == 0.0) ? 1 : 0); }
   Short_t   GetFlag() const                  { return 500; }
   UShort_t  GetNHits() const                 { return (mNHits>0) ? (UShort_t)mNHits : (UShort_t)(-1*mNHits); }
   UShort_t  GetNHitsFit() const              { return (UShort_t)mNHitsFit; }
@@ -104,7 +104,7 @@ class StFemtoTrack : public TObject {
  private:
 
   UShort_t  mId;
-  Char_t    mNHits;           //q * NHitsx
+  Char_t    mNHits;           //q * NHits
   UChar_t   mNHitsPoss;
   UChar_t   mNHitsFit;
   UChar_t   mNHitsDedx;
@@ -129,7 +129,7 @@ class StFemtoTrack : public TObject {
 
   Float_t   mBeta;            // ToF beta
 
-  ClassDef(StFemtoTrack, 1)
+  ClassDef(StFemtoTrack, 2)
 };
 
 #endif

+ 106 - 86
StRoot/StFemtoDstMaker/StHbtFemtoDstReader.cxx

@@ -21,12 +21,39 @@ StHbtFemtoDstReader::StHbtFemtoDstReader(const char *aDir,
   mTChain     = 0;
   mFemtoEvent = 0;
   mEventIndex = 0;
+  mTrig       = 0;
 
   InitRead(mDir, mFileName, mFilter, mMaxFiles);
   mNEvents = (unsigned int)mTChain->GetEntries();
 }
 
 //_________________
+void StHbtFemtoDstReader::SetTrigger(unsigned int trig, char mode) {
+  unsigned int trigColSize = mTrigIdCol.size();
+  for (unsigned int iTrig = 0; iTrig < trigColSize; iTrig++) {
+    if (mTrigIdCol[iTrig] == trig) {
+      unsigned int mask = 1 << iTrig;
+      if (!mode)
+        mTrig = mTrig & (!mask);
+      else
+        mTrig = mTrig | mask;
+      break;
+    }
+  }
+
+  //
+  // Output trigger mask
+  //
+  char bitMask[33] = {0};
+  for (unsigned int iBit = 0; iBit < 32; iBit++) {
+    bitMask[31 - iBit] = '0' + ((mTrig & (1 << iBit)) >> iBit);
+  }
+  std::cout << "StHbtFemtoDstReader[INFO]: Trigger mask = "
+            << mTrig << " bit = "
+            << bitMask << std::endl;
+}
+
+//_________________
 int StHbtFemtoDstReader::FillChain(TChain *aChain, char *aDir,
 				   const char *aFilter, int aMaxFiles) {
 
@@ -125,8 +152,6 @@ int StHbtFemtoDstReader::InitRead(string aDir, string aFileName,
   }
   mTChain->SetBranchAddress("StFemtoEvent", &mFemtoEvent);
   mNEvents = (unsigned int)mTChain->GetEntries();
-  std::cout << "StHbtFemtoDstReader::InitRead [INFO]"
-	    << " NEvents to read: " << mNEvents << std::endl;
   return lNumFiles;
 }
 
@@ -159,6 +184,7 @@ int StHbtFemtoDstReader::GetNEvents() {
 //_________________
 StHbtEvent *StHbtFemtoDstReader::Read() {
 
+
   if(mDebug) {
     std::cout << "StHbtFemtoDstReader[DEBUG]: mEventIndex = "
 	      << mEventIndex << " mNEvents = "
@@ -184,34 +210,28 @@ StHbtEvent *StHbtFemtoDstReader::Read() {
   }
 
   StHbtEvent *mHbtEvent = NULL;
+  mHbtEvent = new StHbtEvent();
 
   if(mFemtoEvent) {
 
-    mHbtEvent = new StHbtEvent();
 
     //Check trigger Id
     if (!mTrigIdCol.empty()) {
-      bool lTrigPass = false;
       unsigned int lTrigId = mFemtoEvent->GetTriggerId();
-      
-      for (unsigned int i = 0; i < mTrigIdCol.size(); i++) {
-	if (lTrigId == mTrigIdCol[i]) {
-	  lTrigPass = true;
-	  break;
-	}
-      }
-      if (!lTrigPass) return 0;
+      if (((lTrigId & mTrig) == 0) || ((lTrigId & (!mTrig)) != 0))
+        return mHbtEvent;
     }
 
+
     float mMagField = mFemtoEvent->GetMagneticField();
     StThreeVectorF v( mFemtoEvent->GetVertexPositionX(),
-		      mFemtoEvent->GetVertexPositionY(),
-		      mFemtoEvent->GetVertexPositionZ() );
+                      mFemtoEvent->GetVertexPositionY(),
+                      mFemtoEvent->GetVertexPositionZ() );
 
     if(mDebug) {
       std::cout << Form("Vertex: x=%4.2f y=%4.2f z=%4.2f ",v.x(),v.y(),v.z()) << std::endl;
     }
-    
+
     mHbtEvent->SetEventNumber( mFemtoEvent->GetEventId() );
     mHbtEvent->SetRunNumber( mFemtoEvent->GetRunId() );
     mHbtEvent->SetZdcAdcEast( mFemtoEvent->GetZDCe() );
@@ -225,85 +245,85 @@ StHbtEvent *StHbtFemtoDstReader::Read() {
     mHbtEvent->SetCent9( mFemtoEvent->GetCent9() );
     mHbtEvent->SetCent16( mFemtoEvent->GetCent16() );
     mHbtEvent->SetRefMultCorr( mFemtoEvent->GetRefMultCorr() );
-    
+
     TClonesArray *mFemtoTrackArr = mFemtoEvent->GetFemtoTracks();
-    
+
     if (mFemtoTrackArr) {
       int mNFemtoTracks = mFemtoTrackArr->GetEntries();
       mTotalTracks += (float)mNFemtoTracks;
       if (mDebug) {
-	std::cout << "StHbtFemtoDstReader[DEBUG]: mNFemtoTracks = "
-		  << mNFemtoTracks << std::endl;
+        std::cout << "StHbtFemtoDstReader[DEBUG]: mNFemtoTracks = "
+            << mNFemtoTracks << std::endl;
       }
-      
+
       for (int i = 0; i < mNFemtoTracks; i++) {
 
-	StFemtoTrack *mFemtoTrack = (StFemtoTrack *)mFemtoTrackArr->At(i);
-	StHbtTrack *mHbtTrack = new StHbtTrack;
-	Short_t charge = mFemtoTrack->GetCharge();
-
-	mHbtTrack->SetTrackId( mFemtoTrack->GetId() );
-	mHbtTrack->SetTrackType(mFemtoTrack->GetType());
-	mHbtTrack->SetCharge( charge );
-	mHbtTrack->SetFlag( mFemtoTrack->GetFlag() );
-	
-	mHbtTrack->SetNHits( mFemtoTrack->GetNHits() );
-	mHbtTrack->SetNHitsFit( mFemtoTrack->GetNHitsFit() );
-	mHbtTrack->SetNHitsPossible( mFemtoTrack->GetNHitsPoss() );
-	mHbtTrack->SetNHitsDedx( mFemtoTrack->GetNHitsDedx() );
-	
-	mHbtTrack->SetNSigmaElectron(mFemtoTrack->GetNSigmaElectron());
-	mHbtTrack->SetNSigmaPion(mFemtoTrack->GetNSigmaPion());
-	mHbtTrack->SetNSigmaKaon(mFemtoTrack->GetNSigmaKaon());
-	mHbtTrack->SetNSigmaProton(mFemtoTrack->GetNSigmaProton());
-	
-	mHbtTrack->SetPidProbElectron(mFemtoTrack->GetPidProbElectron());
-	mHbtTrack->SetPidProbPion(mFemtoTrack->GetPidProbPion());
-	mHbtTrack->SetPidProbKaon(mFemtoTrack->GetPidProbKaon());
-	mHbtTrack->SetPidProbProton(mFemtoTrack->GetPidProbProton());
-	
-	mHbtTrack->SetdEdx( mFemtoTrack->GetDedx() );
-	mHbtTrack->SetChiSquaredXY( mFemtoTrack->GetChi2() );
-	mHbtTrack->SetChiSquaredZ( 1. );
-	
-	mHbtTrack->SetDCAxy( mFemtoTrack->GetDCAxy() );
-	mHbtTrack->SetDCAz( mFemtoTrack->GetDCAz() );
-	mHbtTrack->SetDCAxyGlobal( mFemtoTrack->GetDCAxyGlobal() );
-	mHbtTrack->SetDCAzGlobal( mFemtoTrack->GetDCAzGlobal() );
-
-	StThreeVectorF mGlobMomentum( mFemtoTrack->GetGlobalPx(),
-				      mFemtoTrack->GetGlobalPy(),
-				      mFemtoTrack->GetGlobalPz() );
-	StThreeVectorF mPrimMomentum( mFemtoTrack->GetPx(),
-				      mFemtoTrack->GetPy(),
-				      mFemtoTrack->GetPz() );
-	mHbtTrack->SetP( mPrimMomentum );
-	//mHbtTrack->SetPt( mFemtoTrack->GetPt() );
-
-	if(mDebug) {
-	  std::cout << Form("Track parameters: px=%3.2f py=%3.2f pz=%3.2f q=%1d nhits=%2d",
-			    mPrimMomentum.x(),mPrimMomentum.y(),mPrimMomentum.z(),
-			    charge, mFemtoTrack->GetNHits() ) << std::endl;
-	}
-
-	StThreeVectorF mGlobOrigin( (mFemtoTrack->GetDCAxGlobal() + mFemtoEvent->GetVertexPositionX()),
-				    (mFemtoTrack->GetDCAyGlobal() + mFemtoEvent->GetVertexPositionY()),
-				    (mFemtoTrack->GetDCAzGlobal() + mFemtoEvent->GetVertexPositionZ()) );
-	StThreeVectorF mPrimOrigin( mFemtoEvent->GetVertexPositionX(),
-				    mFemtoEvent->GetVertexPositionY(),
-				    mFemtoEvent->GetVertexPositionZ() );
-	StPhysicalHelixD mGlobHelix( mGlobMomentum, mGlobOrigin, mMagField * kilogauss, charge );
-	StPhysicalHelixD mPrimHelix( mPrimMomentum, mPrimOrigin, mMagField * kilogauss, charge );
-	mHbtTrack->SetHelixGlobal( mGlobHelix );
-	mHbtTrack->SetHelix( mPrimHelix );
-	
-	mHbtTrack->SetTopologyMap( 0, mFemtoTrack->GetMap0() );
-	mHbtTrack->SetTopologyMap( 1, mFemtoTrack->GetMap1() );
-	
-	mHbtTrack->SetTofBeta( mFemtoTrack->GetBeta() );
-	mHbtTrack->SetHiddenInfo( 0 );
-	
-	mHbtEvent->TrackCollection()->push_back( mHbtTrack );
+        StFemtoTrack *mFemtoTrack = (StFemtoTrack *)mFemtoTrackArr->At(i);
+        StHbtTrack *mHbtTrack = new StHbtTrack;
+        Short_t charge = mFemtoTrack->GetCharge();
+
+        mHbtTrack->SetTrackId( mFemtoTrack->GetId() );
+        mHbtTrack->SetTrackType(mFemtoTrack->GetType());
+        mHbtTrack->SetCharge( charge );
+        mHbtTrack->SetFlag( mFemtoTrack->GetFlag() );
+
+        mHbtTrack->SetNHits( mFemtoTrack->GetNHits() );
+        mHbtTrack->SetNHitsFit( mFemtoTrack->GetNHitsFit() );
+        mHbtTrack->SetNHitsPossible( mFemtoTrack->GetNHitsPoss() );
+        mHbtTrack->SetNHitsDedx( mFemtoTrack->GetNHitsDedx() );
+
+        mHbtTrack->SetNSigmaElectron(mFemtoTrack->GetNSigmaElectron());
+        mHbtTrack->SetNSigmaPion(mFemtoTrack->GetNSigmaPion());
+        mHbtTrack->SetNSigmaKaon(mFemtoTrack->GetNSigmaKaon());
+        mHbtTrack->SetNSigmaProton(mFemtoTrack->GetNSigmaProton());
+
+        mHbtTrack->SetPidProbElectron(mFemtoTrack->GetPidProbElectron());
+        mHbtTrack->SetPidProbPion(mFemtoTrack->GetPidProbPion());
+        mHbtTrack->SetPidProbKaon(mFemtoTrack->GetPidProbKaon());
+        mHbtTrack->SetPidProbProton(mFemtoTrack->GetPidProbProton());
+
+        mHbtTrack->SetdEdx( mFemtoTrack->GetDedx() );
+        mHbtTrack->SetChiSquaredXY( mFemtoTrack->GetChi2() );
+        mHbtTrack->SetChiSquaredZ( 1. );
+
+        mHbtTrack->SetDCAxy( mFemtoTrack->GetDCAxy() );
+        mHbtTrack->SetDCAz( mFemtoTrack->GetDCAz() );
+        mHbtTrack->SetDCAxyGlobal( mFemtoTrack->GetDCAxyGlobal() );
+        mHbtTrack->SetDCAzGlobal( mFemtoTrack->GetDCAzGlobal() );
+
+        StThreeVectorF mGlobMomentum( mFemtoTrack->GetGlobalPx(),
+                                     mFemtoTrack->GetGlobalPy(),
+                                     mFemtoTrack->GetGlobalPz() );
+        StThreeVectorF mPrimMomentum( mFemtoTrack->GetPx(),
+                                     mFemtoTrack->GetPy(),
+                                     mFemtoTrack->GetPz() );
+        mHbtTrack->SetP( mPrimMomentum );
+        //mHbtTrack->SetPt( mFemtoTrack->GetPt() );
+
+        if(mDebug) {
+          std::cout << Form("Track parameters: px=%3.2f py=%3.2f pz=%3.2f q=%1d nhits=%2d",
+                            mPrimMomentum.x(),mPrimMomentum.y(),mPrimMomentum.z(),
+                            charge, mFemtoTrack->GetNHits() ) << std::endl;
+        }
+
+        StThreeVectorF mGlobOrigin( (mFemtoTrack->GetDCAxGlobal() + mFemtoEvent->GetVertexPositionX()),
+                                   (mFemtoTrack->GetDCAyGlobal() + mFemtoEvent->GetVertexPositionY()),
+                                   (mFemtoTrack->GetDCAzGlobal() + mFemtoEvent->GetVertexPositionZ()) );
+        StThreeVectorF mPrimOrigin( mFemtoEvent->GetVertexPositionX(),
+                                   mFemtoEvent->GetVertexPositionY(),
+                                   mFemtoEvent->GetVertexPositionZ() );
+        StPhysicalHelixD mGlobHelix( mGlobMomentum, mGlobOrigin, mMagField * kilogauss, charge );
+        StPhysicalHelixD mPrimHelix( mPrimMomentum, mPrimOrigin, mMagField * kilogauss, charge );
+        mHbtTrack->SetHelixGlobal( mGlobHelix );
+        mHbtTrack->SetHelix( mPrimHelix );
+
+        mHbtTrack->SetTopologyMap( 0, mFemtoTrack->GetMap0() );
+        mHbtTrack->SetTopologyMap( 1, mFemtoTrack->GetMap1() );
+
+        mHbtTrack->SetTofBeta( mFemtoTrack->GetBeta() );
+        mHbtTrack->SetHiddenInfo( 0 );
+
+        mHbtEvent->TrackCollection()->push_back( mHbtTrack );
       } //for (int i = 0; i < mNFemtoTracks; i++)
     } //if (mFemtoTrackArr)
   } //if(mFemtoEvent)

+ 3 - 2
StRoot/StFemtoDstMaker/StHbtFemtoDstReader.h

@@ -40,9 +40,9 @@ class StHbtFemtoDstReader : public StHbtEventReader {
   float         mTotalTracks;
 
   std::vector<unsigned int> mTrigIdCol;
+  unsigned int              mTrig;
   
   int mMaxFiles;
-  
 
   int InitRead(string aDir, string aFileName,
 	       string aFilter, int aMaxFiles);
@@ -57,13 +57,14 @@ class StHbtFemtoDstReader : public StHbtEventReader {
 		      const char *aFilter = ".", int aMaxFiles = 1e9, bool aDebug = false);
   StHbtEvent *ReturnHbtEvent();
   void AddTrigId(unsigned int aTrigId);
+  void SetTrigger(unsigned int trig, char mode);
   void PrintTotalTracks();
   void DoDebug(bool aDebug = false);
   StHbtString Report();
   int GetNEvents();
   
 
-  ClassDef(StHbtFemtoDstReader, 2)
+  ClassDef(StHbtFemtoDstReader, 3)
 };
 
 #endif

+ 8 - 8
StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_MomCons.cxx

@@ -54,7 +54,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
   for(int i=0; i<mNumberCFs; i++) {
 
     TitCurrent.str("");
-    TitCurrent << TitNum << "_bin_" << i;
+    TitCurrent << TitNum.str().c_str() << "_bin_" << i;
     mNumerator[i].SetName(TitCurrent.str().c_str());
     mNumerator[i].SetTitle(TitCurrent.str().c_str());
     mNumerator[i].SetBins(nbins,QinvLo,QinvHi);
@@ -62,7 +62,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumerator[i].Sumw2();
     
     TitCurrent.str("");
-    TitCurrent << TitNumPt1xPt2 << "_bin_" << i;
+    TitCurrent << TitNumPt1xPt2.str().c_str() << "_bin_" << i;
     mNumeratorPt1xPt2[i].SetName(TitCurrent.str().c_str());
     mNumeratorPt1xPt2[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorPt1xPt2[i].SetBins(nbins,QinvLo,QinvHi);
@@ -70,7 +70,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumeratorPt1xPt2[i].Sumw2();
 
     TitCurrent.str("");
-    TitCurrent << TitNumPt1SqrSumPt2Sqr << "_bin_" << i;
+    TitCurrent << TitNumPt1SqrSumPt2Sqr.str().c_str() << "_bin_" << i;
     mNumeratorPt1SqrSumPt2Sqr[i].SetName(TitCurrent.str().c_str());
     mNumeratorPt1SqrSumPt2Sqr[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorPt1SqrSumPt2Sqr[i].SetBins(nbins,QinvLo,QinvHi);
@@ -78,7 +78,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumeratorPt1SqrSumPt2Sqr[i].Sumw2();
 
     TitCurrent.str("");
-    TitCurrent << TitNumPz1xPz2 << "_bin_" << i;
+    TitCurrent << TitNumPz1xPz2.str().c_str() << "_bin_" << i;
     mNumeratorPz1xPz2[i].SetName(TitCurrent.str().c_str());
     mNumeratorPz1xPz2[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorPz1xPz2[i].SetBins(nbins,QinvLo,QinvHi);
@@ -86,7 +86,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumeratorPz1xPz2[i].Sumw2();
 
     TitCurrent.str("");
-    TitCurrent << TitNumPz1SqrSumPz2Sqr << "_bin_" << i;
+    TitCurrent << TitNumPz1SqrSumPz2Sqr.str().c_str() << "_bin_" << i;
     mNumeratorPz1SqrSumPz2Sqr[i].SetName(TitCurrent.str().c_str());
     mNumeratorPz1SqrSumPz2Sqr[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorPz1SqrSumPz2Sqr[i].SetBins(nbins,QinvLo,QinvHi);
@@ -94,7 +94,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumeratorPz1SqrSumPz2Sqr[i].Sumw2();
 
     TitCurrent.str("");
-    TitCurrent << TitDen << "_bin_" << i;
+    TitCurrent << TitDen.str().c_str() << "_bin_" << i;
     mDenominator[i].SetName(TitCurrent.str().c_str());
     mDenominator[i].SetTitle(TitCurrent.str().c_str());
     mDenominator[i].SetBins(nbins,QinvLo,QinvHi);
@@ -102,7 +102,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mDenominator[i].Sumw2();
   
     TitCurrent.str("");
-    TitCurrent << TitNumE1xE2 << "_bin_" << i;
+    TitCurrent << TitNumE1xE2.str().c_str() << "_bin_" << i;
     mNumeratorE1xE2[i].SetName(TitCurrent.str().c_str());
     mNumeratorE1xE2[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorE1xE2[i].SetBins(nbins,QinvLo,QinvHi);
@@ -110,7 +110,7 @@ QinvCorrFctnKt_MomCons::QinvCorrFctnKt_MomCons(char* title,
     mNumeratorE1xE2[i].Sumw2();
 
     TitCurrent.str("");
-    TitCurrent << TitNumE1sumE2 << "_bin_" << i;
+    TitCurrent << TitNumE1sumE2.str().c_str() << "_bin_" << i;
     mNumeratorE1sumE2[i].SetName(TitCurrent.str().c_str());
     mNumeratorE1sumE2[i].SetTitle(TitCurrent.str().c_str());
     mNumeratorE1sumE2[i].SetBins(nbins,QinvLo,QinvHi);

+ 1 - 1
StRoot/StHbtMaker/CorrFctn/QinvCorrFctnKt_MomCons.h

@@ -6,7 +6,7 @@
 
 //_________________
 class QinvCorrFctnKt_MomCons : public StHbtCorrFctn {
-
+ public:
   QinvCorrFctnKt_MomCons(char *title,
 			 const int& nbins, const float& QinvLo, const float& QinvHi,
 			 const int& nCFs=20, const float& KtLo=0., const float& KtHi=1.);

+ 1 - 0
StRoot/StHbtMaker/Infrastructure/StHbtManager.cxx

@@ -255,6 +255,7 @@ StHbtEventWriter* StHbtManager::EventWriter( int n ){  // return pointer to n-th
   }
   return *iter;
 }
+
  //____________________________
 int StHbtManager::ProcessEvent(){
 

+ 848 - 0
StRoot/StMuDstQAMaker/StMuDstMinvMaker.cxx

@@ -0,0 +1,848 @@
+#include "StMuDstMinvMaker.h"
+#include "StBTofHeader.h"
+#include "TBranch.h"
+
+ClassImp(StMuDstMinvMaker)
+
+  //
+  // Set maximum file size to 1.9 GB (Root has a 2GB limit)
+  //
+#define MAXFILESIZE 1900000000
+
+
+
+//_________________
+StMuDstMinvMaker::StMuDstMinvMaker(StMuDstMaker *muMaker,
+                                   const Char_t *oFileName) {
+
+  std::cout << "StMuDstMinvMaker::StMuDstMinvMaker - Creating an instance...";
+  mMuDstMaker = muMaker;
+
+  mEventIsGood = false;
+  mNEventsIn = 0;
+  mNEventsPassed = 0;
+
+  mIsGoodTrack = false;
+  mCurrentTrigger = 0;
+
+  mFileName = oFileName;
+  //
+  // Initialize event cut variables
+  //
+  mPrimVtxZ[0] = -70.;
+  mPrimVtxZ[1] = 70.;
+  mPrimVtxR[0] = -1.;
+  mPrimVtxR[1] = 10.;
+  mPrimVtxVpdVzDiff[0] = -10.;
+  mPrimVtxVpdVzDiff[1]= 10.;
+  mPrimVtxXShift = 0.;
+  mPrimVtxYShift = 0.;
+  //
+  // Initialize single-particle cut variables
+  //
+  mTrackP[0] = 0.1;
+  mTrackP[1] = 2.;
+  mTrackDca[0] = 0.;
+  mTrackDca[1] = 5.;
+  mTrackDcaGlobal[0] = 0.;
+  mTrackDcaGlobal[1] = 5.;
+  mTrackNHits[0] = 15;
+  mTrackNHits[1] = 50;
+  mTrackNHitsFit[0] = 15;
+  mTrackNHitsFit[1] = 50;
+
+  mTrackEta[0] = -1.1;
+  mTrackEta[1] = 1.1;
+  mTrackFlag[0] = 0;
+  mTrackFlag[1] = 1000;
+
+  /*
+  //
+  // Initialize TPC and TOF PID
+  //
+  mPionPionNSigma[0] = -3.;   mPionPionNSigma[1] = 3.;
+  mPionKaonNSigma[0] = 1.;    mPionKaonNSigma[1] = -1.;
+  */
+  mTTTPThreshold = 0.7;
+
+  std::cout << "\t[DONE]" << std::endl;
+}
+
+//_________________
+StMuDstMinvMaker::~StMuDstMinvMaker() {
+  /* nothing to do */
+}
+
+//_________________
+Int_t StMuDstMinvMaker::Init() {
+
+  std::cout << "StMuDstMinvMaker::Init - Initializing the maker"
+      << std::endl;
+
+  //
+  // Create histograms and output file
+  //
+  mOutFile = new TFile(mFileName, "RECREATE");
+  //
+  // General event distributions
+  //
+  hR = new TH1F("hR", "R=#sqrt{v_{x}^{2} + v_{y}^{2}};R (cm);#",
+                200, 0., 2.);
+  hVx = new TH1F("hVx", ";v_{x} (cm);#", 200, -1.5, 1.5);
+  hVy = new TH1F("hVy", ";v_{y} (cm);#", 200, -1.5, 1.5);
+  hVz = new TH1F("hVz", ";v_{z} (cm);#", 100, -70., 70.);
+  hVxVsVy = new TH2F("hVxVsVy", ";v_{x} (cm); v_{y} (cm)",
+                     600, -1.5, 1.5,
+                     600, -1.5, 1.5);
+  hVxVsVz = new TH2F("hVxVsVz", ";v_{x} (cm); v_{z} (cm)",
+                     600, -1.5, 1.5,
+                     600, -70., 70.);
+  hVyVsVz = new TH2F("hVyVsVz", ";v_{y} (cm); v_{z} (cm)",
+                     600, -1.5, 1.5,
+                     600, -70., 70.);
+  hNPrimTr = new TH1F("hNPrimTr", ";N_{primary tracks};#", 150, 0, 1500);
+  hNGlobTr = new TH1F("hNGlobTr", ";N_{global tracks};#", 150, 0, 1500);
+  hTofRefMult = new TH1F("hTofRefMult", ";TofRefMult;#", 100, 0, 1000);
+  hTofRefMultVsRefMult = new TH2F("hTofRefMultVsRefMult", ";TofRefMult;RefMult",
+                                  500, 0., 1000,
+                                  500, 0., 1000);
+  hVpd = new TH1F("hVpd", ";v_{z}^{vpd} (cm);#", 100, -70., 70.);
+  hVpdVz = new TH1F("hVpdVz", ";v_{z}^{vpd}-v_{z}", 100, -10., 10.);
+  hNPrimVtx = new TH1F("hNPrimVtx", ";N_{primary vertex};#", 10, 0, 10);
+  //
+  // General track distributions
+  //
+  hP = new TH1F("hP", ";p (GeV/c);#", 200, 0., 5.);
+  hPt = new TH1F("hPt", ";p_{t} (GeV/c);#", 200, 0., 5.);
+  hPx = new TH1F("hPx", ";p_{x} (GeV/c);#", 200, 0., 5.);
+  hPy = new TH1F("hPy", ";p_{y} (GeV/c);#", 200, 0., 5.);
+  hPz = new TH1F("hPz", ";p_{z} (GeV/c);#", 200, 0., 5.);
+  hEta = new TH1F("hEta", ";#eta;#", 200, -2., 2.);
+  hPAccept = new TH1F("hPAccept", "accepted tracks;p (GeV/c);#", 200, 0., 3.);
+  hPtAccept = new TH1F("hPtAccept", "accepted tracks;p_{t} (GeV/c);#", 200, 0., 3.);
+  hPxAccept = new TH1F("hPxAccept", "accepted tracks;p_{x} (GeV/c);#", 200, 0., 3.);
+  hPyAccept = new TH1F("hPyAccept", "accepted tracks;p_{y} (GeV/c);#", 200, 0., 3.);
+  hPzAccept = new TH1F("hPzAccept", "accepted tracks;p_{z} (GeV/c);#", 200, 0., 3.);
+  hEtaAccept = new TH1F("hEtaAccept", "accepted tracks;#eta;#", 200, -2., 2.);
+  hMassSqrVsPt = new TH2F("hMassSqrVsPt", ";p_{t}Q (GeV/c);m^{2} (GeV/c^{2})",
+                          400, -2., 2.,
+                          1000, -0.1, 1.5);
+  hDedxVsPt = new TH2F("hDedxVsPt", ";p_{t}Q (GeV/c);dE/dx (keV/cm)",
+                       400, -2., 2.,
+                       400, 0., 20.);
+  hInvBetaExpVsPt = new TH2F("hInvBetaExpVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}",
+                             400, -2., 2.,
+                             200, 0., 2.);
+  hInvBetaThVsPt = new TH2F("hInvBetaThVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}",
+                            400, -2., 2.,
+                            200, 0., 2.);
+  hTOF = new TH1F("hTOF", "time of flight;t ns;#", 100, 0., 50.);
+  // Pion PID
+  hNSigmaPionPion = new TH1F("hNSigmaPionPion", "#pi;n#sigma(#pi^{#pm});#",
+                             100, -5., 5.);
+  hNSigmaPionKaon = new TH1F("hNSigmaPionKaon", "#pi;n#sigma(#K^{#pm});#",
+                             100, -5., 5.);
+  hNSigmaPionProton = new TH1F("hNSigmaPionProton", "#pi;n#sigma(p);#",
+                               100, -5., 5.);
+  hNSigmaPionPionVsPt = new TH2F("hNSigmaPionPionVsPt",
+                                 "#pi;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
+                                 400, -2., 2.,
+                                 1000, -5., 5.);
+  hNSigmaPionKaonVsPt = new TH2F("hNSigmaPionKaonVsPt",
+                                 "#pi;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
+                                 400, -2., 2.,
+                                 1000, -5., 5.);
+  hNSigmaPionProtonVsPt = new TH2F("hNSigmaPionProtonVsPt",
+                                   "#pi;p_{t}Q (GeV/c);n#sigma(p)",
+                                   400, -2., 2.,
+                                   1000, -5., 5.);
+  hDMassSqrVsPtPionTPC = new TH2F("hDMassSqrVsPtPionTPC",
+                                  "TPC PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})",
+                                  400, -2., 2.,
+                                  500, -.5, .5);
+  hDMassSqrVsPtPionTOF = new TH2F("hDMassSqrVsPtPionTOF",
+                                  "TOF PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})",
+                                  400, -2., 2.,
+                                  500, -.5, .5);
+  hInvBetaExpThPionTNT = new TH2F("hInvBetaExpThPionTNT", "TPC and TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThPionTPC = new TH2F("hInvBetaExpThPionTPC", "only TPC #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThPionTOF = new TH2F("hInvBetaExpThPionTOF", "only TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThPionTTT = new TH2F("hInvBetaExpThPionTTT",
+                                  Form("if p #leq %.2f then TPC, else TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                       mTTTPThreshold),
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThPionINT = new TH2F("hInvBetaExpThPionINT",
+                                  "#pi^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -5., 5.);
+  // Kaon PID
+  hNSigmaKaonPion = new TH1F("hNSigmaKaonPion", "K;n#sigma(#pi^{#pm}};#",
+                             100, -5., 5.);
+  hNSigmaKaonKaon = new TH1F("hNSigmaKaonKaon", "K;n#sigma(#K^{#pm}};#",
+                             100, -5., 5.);
+  hNSigmaKaonProton = new TH1F("hNSigmaKaonProton", "K;n#sigma(p);#",
+                               100, -5., 5.);
+  hNSigmaKaonPionVsPt = new TH2F("hNSigmaKaonPionVsPt",
+                                 "K;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
+                                 400, -2., 2.,
+                                 1000, -5., 5.);
+  hNSigmaKaonKaonVsPt = new TH2F("hNSigmaKaonKaonVsPt",
+                                 "K;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
+                                 400, -2., 2.,
+                                 1000, -5., 5.);
+  hNSigmaKaonProtonVsPt = new TH2F("hNSigmaKaonProtonVsPt",
+                                   "K;p_{t}Q (GeV/c);n#sigma(p)",
+                                   400, -2., 2.,
+                                   1000, -5., 5.);
+  hDMassSqrVsPtKaonTPC = new TH2F("hDMassSqrVsPtKaonTPC",
+                                  "TPC PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})",
+                                  400, -2., 2.,
+                                  500, -.5, .5);
+  hDMassSqrVsPtKaonTOF = new TH2F("hDMassSqrVsPtKaonTOF",
+                                  "TOF PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})",
+                                  400, -2., 2.,
+                                  500, -.5, .5);
+  hInvBetaExpThKaonTNT = new TH2F("hInvBetaExpThKaonTNT", "TPC and TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThKaonTPC = new TH2F("hInvBetaExpThKaonTPC", "only TPC K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThKaonTOF = new TH2F("hInvBetaExpThKaonTOF", "only TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThKaonTTT = new TH2F("hInvBetaExpThKaonTTT",
+                                  Form("if p #leq %.2f then TPC, else TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                       mTTTPThreshold),
+                                  400, -2., 2.,
+                                  500, -0.5, 0.5);
+  hInvBetaExpThKaonINT = new TH2F("hInvBetaExpThKaonINT",
+                                  "K^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                  400, -2., 2.,
+                                  500, -5., 5.);
+  // Proton PID
+  hNSigmaProtonPion = new TH1F("hNSigmaProtonPion", "p;n#sigma(#pi^{#pm}};#",
+                               100, -5., 5.);
+  hNSigmaProtonKaon = new TH1F("hNSigmaProtonKaon", "p;n#sigma(#K^{#pm}};#",
+                               100, -5., 5.);
+  hNSigmaProtonProton = new TH1F("hNSigmaProtonProton", "p;n#sigma(p);#",
+                                 100, -5., 5.);
+  hNSigmaProtonPionVsPt = new TH2F("hNSigmaProtonPionVsPt",
+                                   "p;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
+                                   400, -2., 2.,
+                                   1000, -5., 5.);
+  hNSigmaProtonKaonVsPt = new TH2F("hNSigmaProtonKaonVsPt",
+                                   "p;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
+                                   400, -2., 2.,
+                                   1000, -5., 5.);
+  hNSigmaProtonProtonVsPt = new TH2F("hNSigmaProtonProtonVsPt",
+                                     "p;p_{t}Q (GeV/c);n#sigma(p)",
+                                     400, -2., 2.,
+                                     1000, -5., 5.);
+  hDMassSqrVsPtProtonTPC = new TH2F("hDMassSqrVsPtProtonTPC",
+                                    "TPC PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})",
+                                    400, -2., 2.,
+                                    500, -.5, .5);
+  hDMassSqrVsPtProtonTOF = new TH2F("hDMassSqrVsPtProtonTOF",
+                                    "TOF PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})",
+                                    400, -2., 2.,
+                                    500, -.5, .5);
+  hInvBetaExpThProtonTNT = new TH2F("hInvBetaExpThProtonTNT", "TPC and TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                    400, -2., 2.,
+                                    500, -0.5, 0.5);
+  hInvBetaExpThProtonTPC = new TH2F("hInvBetaExpThProtonTPC", "only TPC p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                    400, -2., 2.,
+                                    500, -0.5, 0.5);
+  hInvBetaExpThProtonTOF = new TH2F("hInvBetaExpThProtonTOF", "only TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                    400, -2., 2.,
+                                    500, -0.5, 0.5);
+  hInvBetaExpThProtonTTT = new TH2F("hInvBetaExpThProtonTTT",
+                                    Form("if p #leq %.2f then TPC, else TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                         mTTTPThreshold),
+                                    400, -2., 2.,
+                                    500, -0.5, 0.5);
+  hInvBetaExpThProtonINT = new TH2F("hInvBetaExpThProtonINT",
+                                    "p;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
+                                    400, -2., 2.,
+                                    500, -5., 5.);
+  //
+  // Any charge
+  //
+  hRefMult = new TH1F("hRefMult", ";RefMult;#", 100, 0., 1000.);
+  hRefMultAccept = new TH1F("hRefMultAccept", ";RefMultAccept;#", 100, 0., 1000.);
+
+  std::cout << "StMuDstMinvMaker::Init - Initialization has been finished"
+      << std::endl;
+  return StMaker::Init();
+}
+
+//________________
+void StMuDstMinvMaker::Clear(Option_t *option) {
+  StMaker::Clear();
+}
+
+//________________
+Int_t StMuDstMinvMaker::Make() {
+
+  mNEventsIn++;
+
+  mMuDst = NULL;
+  mMuEvent = NULL;
+
+  //MuDstMaker initialization
+  mMuDstMaker = (StMuDstMaker*)GetMaker("MuDst");
+  if(!mMuDstMaker) {
+    LOG_ERROR << "StMuDstMinvMaker::Make [ERROR] - Cannot find StMuDstMaker"
+        << std::endl;
+    return kStOk;
+  }
+
+  //Obtaining MuDst
+  mMuDst = (StMuDst*)GetInputDS("MuDst");
+  if(!mMuDst) {
+    gMessMgr->Warning() << "StMuDstMinvMaker::Make [WARNING] - No MuDst has been found" 
+        << endm;
+    return kStOk;
+  }
+
+  //Obtaining MuEvent
+  mMuEvent = (StMuEvent*)mMuDst->event();
+  if(!AcceptTrigger(mMuEvent)) { //If trigger is not found
+    return kStOk;
+  }
+
+  //Multiplicity cannot be negative
+  unsigned short refMult = mMuEvent->refMult();
+  unsigned short refMultPos = 0;
+  unsigned short refMultNeg = 0;
+
+  if(refMult < 0) {
+    return kStOk;
+  }
+
+  hRefMult->Fill(refMult);
+
+  int mNVertices = mMuDst->numberOfPrimaryVertices();
+  int mNPrimTracks = mMuDst->numberOfPrimaryTracks();
+  int mNGlobTracks = mMuDst->numberOfGlobalTracks();
+
+  //Some initializations of local variables
+  mPrimVertex = NULL;
+  StThreeVectorF mVertPosition;
+  Float_t mVpdVz = 0.;
+  Int_t   mPrimVertIndex = -999;
+  Float_t mRanking = -999.;
+
+  //Clean index vectors
+  CleanVariables();
+
+  //Vertex loop
+  unsigned short refMultAccept = 0;
+  unsigned short tofRefMult = 0;
+  unsigned short refMultAcceptPos = 0;
+  unsigned short refMultAcceptNeg = 0;
+  for(Int_t iVert=0; iVert<mNVertices; iVert++) {
+
+    mEventIsGood = false;
+
+    hNPrimVtx->Fill(mNVertices);
+
+    if(iVert!=0) continue; // Not first primary vertex does not contain tracks with fast detectors (TOF)
+
+    mPrimVertex = mMuDst->primaryVertex(iVert);
+    //Positive ranking only
+    if(mPrimVertex->ranking() <= 0) continue;
+    //Not (0,0,0) position of the primary vertex
+    if(mPrimVertex->position().x() == 0 &&
+       mPrimVertex->position().y() == 0 &&
+       mPrimVertex->position().z() == 0) continue;
+    //Reasonable amount of tracks
+    if(mNPrimTracks < 0 || mNPrimTracks > 10000) continue;
+
+    mPrimVertIndex = iVert;
+    mRanking = mPrimVertex->ranking();
+    mVertPosition = mPrimVertex->position();
+    mVpdVz = mMuDst->btofHeader()->vpdVz();
+
+    if( !AcceptPrimVtx(mVertPosition, mVpdVz) ) continue;
+
+    mEventIsGood = true;
+
+    hNPrimTr->Fill(mNPrimTracks);
+    hNGlobTr->Fill(mNGlobTracks);
+    //
+    //Loop over primary tracks
+    //
+    for(int iTrk = 0; iTrk < mNPrimTracks; iTrk++) {
+      mPrimTrack = mMuDst->primaryTracks(iTrk);
+      mGlobTrack = mMuDst->globalTracks(mPrimTrack->index2Global());
+      short charge = mPrimTrack->charge() > 0 ? 1 : -1;
+      float eta = mPrimTrack->eta();
+      float pt  = mPrimTrack->pt();
+      float p   = mPrimTrack->p().mag();
+      float p2  = mPrimTrack->p().mag2();
+      float px  = mPrimTrack->p().x();
+      float py  = mPrimTrack->p().y();
+      float pz  = mPrimTrack->p().z();
+
+      charge > 0 ? refMultPos++ : refMultNeg++;
+
+      hP->Fill(p);
+      hPt->Fill(pt);
+      hPx->Fill(px);
+      hPy->Fill(py);
+      hPz->Fill(pz);
+      hEta->Fill(eta);
+
+      if( !AcceptTrack(mPrimTrack, iVert) ) continue;
+
+      hPAccept->Fill(p);
+      hPtAccept->Fill(pt);
+      hPxAccept->Fill(px);
+      hPyAccept->Fill(py);
+      hPzAccept->Fill(pz);
+      hEtaAccept->Fill(eta);
+
+      refMultAccept++;
+      charge > 0 ? refMultAcceptPos++ : refMultAcceptNeg++;
+
+      float nSigmaPion = mPrimTrack->nSigmaPion();
+      float nSigmaKaon = mPrimTrack->nSigmaKaon();
+      float nSigmaProton = mPrimTrack->nSigmaProton();
+      float dedx = mPrimTrack->dEdx();
+      bool  tofTrack = IsTofTrack(mGlobTrack);
+      float betaExp;
+      float massSqr;
+      bool pionTPC   = ((nSigmaPion   >= mPionPionNSigma[0]     || nSigmaPion   <= mPionPionNSigma[1])      && 
+                        (nSigmaKaon   <  mPionKaonNSigma[0]     || nSigmaKaon   >  mPionKaonNSigma[1])      &&
+                        (nSigmaProton <  mPionProtonNSigma[0]   || nSigmaProton >  mPionProtonNSigma[1]));
+      bool kaonTPC   = ((nSigmaPion   <  mKaonPionNSigma[0]     || nSigmaPion   >  mKaonPionNSigma[1])      &&
+                        (nSigmaKaon   >= mKaonKaonNSigma[0]     || nSigmaKaon   <= mKaonKaonNSigma[1])      &&
+                        (nSigmaProton <  mKaonProtonNSigma[0]   || nSigmaProton >  mKaonProtonNSigma[1]));
+      bool protonTPC = ((nSigmaPion   <  mProtonPionNSigma[0]   || nSigmaPion   >  mProtonPionNSigma[1])    &&
+                        (nSigmaKaon   <  mProtonKaonNSigma[0]   || nSigmaKaon   >  mProtonKaonNSigma[1])    &&
+                        (nSigmaProton >= mProtonProtonNSigma[0] || nSigmaProton <= mProtonProtonNSigma[1]));
+      bool pionTOF   = false;
+      bool kaonTOF   = false;
+      bool protonTOF = false;
+      //
+      // Pion TPC
+      //
+      if (pionTPC) {
+        hNSigmaPionPion->Fill(nSigmaPion);
+        hNSigmaPionKaon->Fill(nSigmaKaon);
+        hNSigmaPionProton->Fill(nSigmaProton);
+        hNSigmaPionPionVsPt->Fill(pt*charge, nSigmaPion);
+        hNSigmaPionKaonVsPt->Fill(pt*charge, nSigmaKaon);
+        hNSigmaPionProtonVsPt->Fill(pt*charge, nSigmaProton);
+      }
+      //
+      // Kaon TPC
+      //
+      if (kaonTPC) {
+        hNSigmaKaonPion->Fill(nSigmaPion);
+        hNSigmaKaonKaon->Fill(nSigmaKaon);
+        hNSigmaKaonProton->Fill(nSigmaProton);
+        hNSigmaKaonPionVsPt->Fill(pt*charge, nSigmaPion);
+        hNSigmaKaonKaonVsPt->Fill(pt*charge, nSigmaKaon);
+        hNSigmaKaonProtonVsPt->Fill(pt*charge, nSigmaProton);
+      }
+      //
+      // Proton TPC
+      //
+      if (protonTPC) {
+        hNSigmaProtonPion->Fill(nSigmaPion);
+        hNSigmaProtonKaon->Fill(nSigmaKaon);
+        hNSigmaProtonProton->Fill(nSigmaProton);
+        hNSigmaProtonPionVsPt->Fill(pt*charge, nSigmaPion);
+        hNSigmaProtonKaonVsPt->Fill(pt*charge, nSigmaKaon);
+        hNSigmaProtonProtonVsPt->Fill(pt*charge, nSigmaProton);
+      }
+
+      hDedxVsPt->Fill(pt*charge, dedx*1e6);
+
+      if (tofTrack) {
+        hTOF->Fill(mGlobTrack->btofPidTraits().timeOfFlight());
+
+        float betaThPion = sqrt(p2/(PION_MASS*PION_MASS + p2));
+        float betaThKaon = sqrt(p2/(KAON_MASS*KAON_MASS + p2));
+        float betaThProton = sqrt(p2/(PROTON_MASS*PROTON_MASS + p2));
+        tofRefMult++;
+        betaExp = mGlobTrack->btofPidTraits().beta();
+        massSqr = p2*(1./(betaExp*betaExp) - 1.);
+        hMassSqrVsPt->Fill(pt*charge, massSqr);
+        pionTOF   = massSqr >= mPionMass[0]   && massSqr <= mPionMass[1];
+        kaonTOF   = massSqr >= mKaonMass[0]   && massSqr <= mKaonMass[1];
+        protonTOF = massSqr >= mProtonMass[0] && massSqr <= mProtonMass[1];
+
+        float dMpion   = massSqr - PION_MASS*PION_MASS;
+        float dMkaon   = massSqr - KAON_MASS*KAON_MASS;
+        float dMproton = massSqr - PROTON_MASS*PROTON_MASS;
+        hInvBetaExpVsPt->Fill(pt*charge, 1./betaExp);
+        //
+        // TPC PID
+        //
+        if (pionTPC) {
+          hInvBetaExpThPionTPC->Fill(pt*charge, 1./betaExp - 1./betaThPion);
+          hDMassSqrVsPtPionTPC->Fill(pt*charge, dMpion);
+        }
+        else if (kaonTPC) {
+          hInvBetaExpThKaonTPC->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
+          hDMassSqrVsPtKaonTPC->Fill(pt, dMkaon);
+        }
+        else if (protonTPC) {
+          hInvBetaExpThProtonTPC->Fill(pt*charge, 1./betaExp - 1./betaThProton);
+          hDMassSqrVsPtProtonTPC->Fill(pt, dMproton);
+        }
+        //
+        // TOF PID
+        //