Nikita Ermakov 2 years ago
parent
commit
cba4c065f0

BIN
StRoot/StFemtoDstMaker.tar.gz


+ 19 - 6
StRoot/StFemtoDstMakerSph/StFemtoDstQAMaker.cxx

@@ -37,9 +37,9 @@ Int_t StFemtoDstQAMaker::Init() {
     Int_t mRefMultWeightBins = 200;
     Double_t mRefMultWeightLo = 0.;
     Double_t mRefMultWeightHi = 20.;
-    Int_t mRefMultBins = 700;
-    Double_t mRefMultLo = -0.5;
-    Double_t mRefMultHi = 699.5;
+    Int_t mRefMultBins = 350;
+    Double_t mRefMultLo = 0.;
+    Double_t mRefMultHi = 350.;
     Int_t mPrimVtxZBins = 400;
     Double_t mPrimVtxZLo = -200.;
     Double_t mPrimVtxZHi = 200.;
@@ -96,6 +96,12 @@ Int_t StFemtoDstQAMaker::Init() {
                             mRefMultBins, mRefMultLo, mRefMultHi);
     hTofMult = new TH1F("hTofMult", ";TofMult;#",
                         mRefMultBins, mRefMultLo, mRefMultHi);
+    hTofRefDiff = new TH1F("hTofRefDiff", ";TofMult - RefMult;#",
+                        2*mRefMultBins, -mRefMultHi, mRefMultHi);
+    hTofRef2Diff = new TH1F("hTofRef2Diff", ";TofMult - RefMult2;#",
+                        2*mRefMultBins, -mRefMultHi, mRefMultHi);
+    hRefRef2Diff = new TH1F("hRefRef2Diff", ";RefMult - RefMult2;#",
+                        2*mRefMultBins, -mRefMultHi, mRefMultHi);
     hTofVsRef = new TH2F("hTofVsRef", ";TofMult;RefMult",
                          mRefMultBins, mRefMultLo, mRefMultHi,
                          mRefMultBins, mRefMultLo, mRefMultHi);
@@ -115,14 +121,16 @@ Int_t StFemtoDstQAMaker::Init() {
                              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]);
+    hRefMultVsRunNumber = new TProfile("hRefMultVsRunNumber","refMult vs. Run Number;Run Number;<refMult>",
+                                       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]);
+    hRefMult2VsRunNumber = new TProfile("hRefMult2VsRunNumber","RefMult2 vs. Run Number;Run Number;<RefMult2>",
+                                       mRunIdBins, mRunIdVal[0], mRunIdVal[1]);
     hTriggerId = new TH1I("hTriggerId", ";trigger bit;#", 32, 0, 32);
     hSphericity = new TH1F("hSphericity", ";sphericity;#", 120, -0.1, 1.1);
     // Track histograms
@@ -365,6 +373,9 @@ Int_t StFemtoDstQAMaker::Make() {
         hRefMultWeight->Fill(femtoEvent->GetRefMultCorrWeight());
         hRefMult->Fill(femtoEvent->GetRefMult());
         hRefMult2->Fill(femtoEvent->GetRefMult2());
+        hTofRefDiff->Fill((int)femtoEvent->GetNumberOfBTofHit() - (int)femtoEvent->GetRefMult());
+        hTofRef2Diff->Fill((int)femtoEvent->GetNumberOfBTofHit() - (int)femtoEvent->GetRefMult2());
+        hRefRef2Diff->Fill((int)femtoEvent->GetRefMult() - (int)femtoEvent->GetRefMult2());
         hRefMultCorr->Fill(femtoEvent->GetRefMultCorr());
         hTofMult->Fill(femtoEvent->GetNumberOfBTofHit());
         hTofVsRef->Fill(femtoEvent->GetNumberOfBTofHit(),
@@ -378,13 +389,15 @@ Int_t StFemtoDstQAMaker::Make() {
         hPrimVertXvsY->Fill(femtoEvent->GetVertexPositionX(),
                             femtoEvent->GetVertexPositionY());
         hPrimVertVpdVzDiff->Fill(femtoEvent->GetVertexPositionZ() - femtoEvent->GetVpdVz());
-        hRefMultVsRunNumer->Fill(femtoEvent->GetRunId(),femtoEvent->GetRefMult());
+        hRefMultVsRunNumber->Fill(femtoEvent->GetRunId(),femtoEvent->GetRefMult());
         hPrimTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),
                                        femtoEvent->GetNumberOfPrimaryTracks());
         hGlobTrackNumVsRunNumber->Fill(femtoEvent->GetRunId(),
                                        femtoEvent->GetNumberOfGlobalTracks());
         hTofMultVsRunNumber->Fill(femtoEvent->GetRunId(),
                                   femtoEvent->GetNumberOfBTofHit());
+        hRefMult2VsRunNumber->Fill(femtoEvent->GetRunId(),
+                                   femtoEvent->GetRefMult2());
         hSphericity->Fill(femtoEvent->GetSphericity());
 
         TClonesArray *lTracks = femtoEvent->GetFemtoTracks();

+ 5 - 1
StRoot/StFemtoDstMakerSph/StFemtoDstQAMaker.h

@@ -53,6 +53,9 @@ class StFemtoDstQAMaker : public StMaker {
   TH1F *hRefMult2;
   TH1F *hRefMultWeight;
   TH1F *hTofMult;
+  TH1F *hTofRefDiff;
+  TH1F *hTofRef2Diff;
+  TH1F *hRefRef2Diff;
   TH2F *hTofVsRef;
   TH2F *hTofVsRef2;
   TH2F *hRefVsRef2;
@@ -62,10 +65,11 @@ class StFemtoDstQAMaker : public StMaker {
   TH1F *hPrimVertVpdVzDiff;
   TH1I *hTriggerId;
   //Here is the way of how to look on the quantities over the run
-  TProfile *hRefMultVsRunNumer;
+  TProfile *hRefMultVsRunNumber;
   TProfile *hPrimTrackNumVsRunNumber;
   TProfile *hGlobTrackNumVsRunNumber;
   TProfile *hTofMultVsRunNumber;
+  TProfile *hRefMult2VsRunNumber;
   //Track histograms
   TH2F *hNSigmaPionVsMom;
   TH2F *hNSigmaKaonVsMom;

+ 7 - 5
StRoot/StFemtoDstMakerSph/StFemtoEvent.h

@@ -11,6 +11,8 @@
 #include "StFemtoTrack.h"
 #include "TClonesArray.h"
 #include "TVector3.h"
+#include "math.h"
+
 
 //_________________
 class StFemtoEvent : public TObject {
@@ -43,7 +45,7 @@ class StFemtoEvent : public TObject {
   void SetVpdVz(Float_t vz)           { mVpdVz = vz; } 
   void SetTriggerId(UInt_t id)        { mTriggerId = id; }
   void SetPrimaryVertexRanking(Float_t);
-  void SetSphericity(Float_t sph)     { mSphericity = sph; }
+  void SetSphericity(Float_t sph)     { mSphericity = (char)round(sph*100.); }
   StFemtoTrack* AddFemtoTrack();
   void Clear(Option_t* option = "C");
 
@@ -77,13 +79,13 @@ class StFemtoEvent : public TObject {
   Float_t  GetPrimaryVertexRanking() const {return (Float_t)mPrimaryVertexRanking; }
   TClonesArray *GetFemtoTracks() const { return mFemtoTracks; }
   Int_t    GetNFemtoTracks() const { return (Int_t)mNFemtoTracks; }
-  Float_t  GetSphericity() const { return mSphericity; }
+  Float_t  GetSphericity() const { return ((float)mSphericity)/100.; }
 
  private:
-  
+
   Int_t    mEventId;
   Int_t    mRunId;
-  Float_t  mSphericity;          //transverse sphericity
+  Char_t   mSphericity;          //transverse sphericity
   Bool_t   mCollisionType;       //0-pp, 1-AA(from those where RefMultCorr is defined)
   UShort_t mRefMult;             //mRefMult
   UShort_t mRefMultCorr;         //mRefMultCorr * 10
@@ -106,7 +108,7 @@ class StFemtoEvent : public TObject {
   UShort_t mNFemtoTracks;
   TClonesArray *mFemtoTracks;
 
-  ClassDef(StFemtoEvent, 100)
+  ClassDef(StFemtoEvent, 101)
 };
 
 #endif

+ 31 - 10
StRoot/StSphericity/StSphericity.cxx

@@ -17,7 +17,7 @@ StSphericity::StSphericity(StMuDstMaker *muMaker,
     matrix       = new TMatrixTSym<double>(2);
     matrix3d     = new TMatrixTSym<double>(3);
     sphLoCut     = 0.1;
-    sphHiCut     = 0.8;
+    sphHiCut     = 0.5;
 }
 
 
@@ -36,34 +36,35 @@ Int_t StSphericity::Init()
     //
     // Create histograms
     //
+    hRefMult = new TH1I("hRefMult", ";RefMult;events", 100, 0, 100);
     // no cut to sphericity
     hEta = new TH1F("hEta", ";#eta;#", 100, -3., 3.);
     hPhi = new TH1F("hPhi", ";#phi [rad];#", 300, -7., 7.);
     hDEtaDPhi = new TH2F("hDEtaDPhi", ";#Delta#eta;#Delta#phi [rad]",
                          200, -3., 3., /* eta */
-                         400, -10., 10.  /* phi */);
+                         400, -1.5, 6. /* phi */);
     hEtaPhi = new TH2F("hEtaPhi", ";#eta;#phi [rad]", 
                        200, -3., 3., /* eta */
-                       400, -10., 10.  /* phi */);
+                       400, -1.5, 6. /* phi */);
     // with sphericity cut
     hDEtaDPhi_lo_sphManual = new TH2F("hDEtaDPhi_lo_sphManual", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sphManual = new TH2F("hDEtaDPhi_hi_sphManual", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_lo_sphAuto = new TH2F("hDEtaDPhi_lo_sphAuto", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sphAuto = new TH2F("hDEtaDPhi_hi_sphAuto", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_lo_sph3d = new TH2F("hDEtaDPhi_lo_sph3d", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sph3d = new TH2F("hDEtaDPhi_hi_sph3d", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     // transverse and 3d sphericity histograms
     hSphManual = new TH1F("hSphManual", ";Sphericity;#",
                           300, -0.1, 1.1);
@@ -231,10 +232,15 @@ Int_t StSphericity::Make()
         matrix->Zero();
         matrix3d->Zero();
 
+        hRefMult->Fill(refMult);
+
         for (int i = 0; i < nPrimTracks; i++)
         {
             mPrimTrack = mMuDst->primaryTracks(i);
-            if (mPrimTrack->vertexIndex() != 0)   return kStOk;
+            //
+            // Primary track cut's
+            //
+            if (!AcceptTrack(mPrimTrack))   continue;
 
             float pt = mPrimTrack->pt();
             const StThreeVectorF &p = mPrimTrack->p();
@@ -339,6 +345,7 @@ Int_t StSphericity::Make()
         {
             if (i == numOfLeadPart)   continue;
             StMuTrack *lTrack = mMuDst->primaryTracks(i);
+            if (!AcceptTrack(lTrack))   continue;
             if (lTrack->p().mag() > 0.5)
             {
                 float phi = lTrack->phi();
@@ -370,3 +377,17 @@ Int_t StSphericity::Make()
     }
     return kStOk;
 }
+
+
+
+bool StSphericity::AcceptTrack(StMuTrack *track)
+{
+    return (track->index2Global() >= 0 &&
+            track->vertexIndex() == 0 &&
+            track->type() == 1 &&
+            track->nHits() >= 15 &&
+            track->nHitsFit() >= 15 &&
+            track->flag() >= 300 && track->flag() < 1000 &&
+            track->dca(0).perp() >= 0. && track->dca(0).perp() <= 3. &&
+            track->dcaGlobal(0).perp() >= 0. && track->dcaGlobal(0).perp() <= 3.);
+}

+ 3 - 1
StRoot/StSphericity/StSphericity.h

@@ -57,7 +57,7 @@ class StSphericity : public StMaker
      void SetSphericityLowCut(float val) { sphLoCut = val; }
      void SetSphericityHighCut(float val) { sphHiCut = val; }
      void SetPtCut(float pt) { ptCut = pt; }
-
+     bool AcceptTrack(StMuTrack *track);
 
     private:
      TFile        *mOutFile;
@@ -77,6 +77,8 @@ class StSphericity : public StMaker
 
      std::vector<unsigned int> mTrigIdCol;
 
+     TH1I *hRefMult;
+
      TH1F *hEta;
      TH1F *hPhi;
      TH2F *hEtaPhi;

+ 7 - 9
StRoot/StSphericityStarGen/StSphericity.cxx

@@ -51,29 +51,27 @@ Int_t StSphericity::Init()
     hPhi = new TH1F("hPhi", ";#phi [rad];#", 300, -7., 7.);
     hDEtaDPhi = new TH2F("hDEtaDPhi", ";#Delta#eta;#Delta#phi [rad]",
                          200, -3., 3., /* eta */
-                         400, -10., 10.  /* phi */);
+                         400, -1.5, 6. /* phi */);
     hEtaPhi = new TH2F("hEtaPhi", ";#eta;#phi [rad]", 
                        200, -3., 3., /* eta */
-                       400, -10., 10.  /* phi */);
+                       400, -1.5, 6. /* phi */);
     // with sphericity cut
     hDEtaDPhi_lo_sphManual = new TH2F("hDEtaDPhi_lo_sphManual", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sphManual = new TH2F("hDEtaDPhi_hi_sphManual", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_lo_sphAuto = new TH2F("hDEtaDPhi_lo_sphAuto", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sphAuto = new TH2F("hDEtaDPhi_hi_sphAuto", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_lo_sph3d = new TH2F("hDEtaDPhi_lo_sph3d", ";#Delta#eta;#Delta#phi [rad]",
                                       200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
+                                      400, -1.5, 6.  /* phi */);
     hDEtaDPhi_hi_sph3d = new TH2F("hDEtaDPhi_hi_sph3d", ";#Delta#eta;#Delta#phi [rad]",
-                                      200, -3., 3., /* eta */
-                                      400, -10., 10.  /* phi */);
     // transverse and 3d sphericity histograms
     hSphManual = new TH1F("hSphManual", ";Sphericity;#",
                           300, -0.1, 1.1);

+ 519 - 0
macros/HBT/pp/nomult/expKaonHBT_pp200_y2012_1D.C

@@ -0,0 +1,519 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+//
+// To define cut monitors
+// uncomments line below
+//
+#define MONITORS
+
+//
+// Constants for analyses
+//
+
+#define CHARGE_BINS 2
+#define SPH_BINS 10
+
+//
+//  Event information
+//
+const int mPrimVertZbins = 6;                // number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};    // redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                // distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};  // VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const bool mUseZdcCorrection = true;         // valid for AuAu collisions. True for 200GeV only!
+
+//
+// Track information
+//
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1;                     // primary track
+const int mPythiaKaonCode = 321;              // pythia particle code
+const int mPythiaPionCode = 211;              // pythia particle code
+
+const int mTrackNHits[2] = {10, 50};          // number of track hits
+const double mTrackNHitsFit[2] = {5, 50};    // number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};         // DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.;
+const double mTrackMom[2] = {0.12, 1.55};        // general momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.05, 1.55};  // general transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -0.5, 0.5};
+
+const short mDetectorSelection = 4;  // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  // track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2] = {-1e9, 1e9};
+const double mnSigmaKaon[2] = {-2., 2.};
+const double mnSigmaProton[2] = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.45};   // track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; // based on p+p data
+const double mKaonTofMsqr[2] = {0.2, 0.32};    // based on p+p data
+
+const double mTrackTnTMom[2] = {0.15, 1.55};
+const double mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-1e9, 1e9};
+const double mTnTNSigmaKaon[2] = {-2., 2.};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.5; // TPC&&TOF (P>Pthr)||TPC(P<Pthr)
+
+//
+// Pair information
+//
+const double mFmrCut[2] = {-1.1, 0.1};         // not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6};  // splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {-1e9., 1e9};        // pair transverse momentum cut
+const double mPairQinv[2] = {-1e9., 1e9};      // pair relative momentum
+const double mAveSep[2] = {5., 5000.};         // average separation of tracks within TPC
+
+//
+// Histograms values
+//
+const int mPairKtBins = 28;                // bin width 
+const double mPairKtVal[2] = {0.1, 1.5};   // 50 MeV/c 
+const int mQinvBins = 150;                 // 20 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 3.0};      // redefine for each analysis
+
+const Char_t *defaultInFile = "test.list";
+const Char_t *defaultOutFile = "oExpKaon_pp200_1D.root";
+
+
+
+void expKaonHBT_pp200_y2012_1D(const Char_t *inFileList = defaultInFile,
+                               const Char_t *oFileName  = defaultOutFile)
+{
+    std::cout << "*****************************" << std::endl
+              << "*    expKaonHBT_pp200_1D    *" << std::endl
+              << "*          Start            *" << std::endl
+              << "*****************************" << std::endl << std::endl;
+
+    //
+    // Load libraries
+    //
+    std::cout << "Loading libraries..." << std::endl;
+    gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+    loadSharedLibraries();
+    gSystem->Load("libMinuit");
+    gSystem->Load("StRefMultCorr");
+    gSystem->Load("StFlowMaker");   // should be placed before HBT
+    gSystem->Load("StHbtMaker");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("StPicoDstMakerRun12");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libStDb_Tables.so");
+    gSystem->Load("libgen_Tables.so");
+    gSystem->Load("libgeometry_Tables.so");
+    gSystem->Load("libsim_Tables.so");
+    gSystem->Load("libStarMagField.so");
+    gSystem->Load("libSt_db_Maker.so");
+    gSystem->Load("libSt_g2t.so");
+    gSystem->Load("libSt_geant_Maker.so");
+    gSystem->Load("StarGeneratorUtil.so");
+    gSystem->Load("StarGeneratorEvent.so");
+    gSystem->Load("StarGeneratorBase.so");
+    gSystem->Load("StFemtoDstMakerSph");
+
+    std::cout << "Libraries have been successfully loaded" << std::endl;
+
+    //
+    // Create chain
+    //
+    StChain* mChain = new StChain("StChain");
+    mChain->SetDebug(0);
+    StMuDebug::setLevel(0);
+
+    //
+    // Set StHbtMaker and StHbtManager
+    //
+    StHbtMaker *mHbtMaker = new StHbtMaker();
+    mHbtMaker->SetDebug(0);
+    StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+    //
+    // Set StHbtReader
+    //
+    StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+    femtoReader->AddTrigId(370011); // vpdmb-nobsmd
+    femtoReader->AddTrigId(370001); // vpdmb
+    femtoReader->AddTrigId(370022); // bbcmb
+    femtoReader->AddTrigId(370501); // bht0-vpdmb
+    femtoReader->AddTrigId(370511); // bht1-vpdmb
+    femtoReader->AddTrigId(370531); // bht2
+    femtoReader->AddTrigId(370542); // bht0-bbcmb-tof0
+    femtoReader->AddTrigId(370546); // bht1-bbcmb-tof0
+    femtoReader->AddTrigId(370522); // bht2-bbcmb
+    femtoReader->AddTrigId(370601); // jp0
+    femtoReader->AddTrigId(370611); // jp1
+    femtoReader->AddTrigId(370621); // jp2
+    femtoReader->AddTrigId(370641); // ajp
+    femtoReader->AddTrigId(370301); // bbcmb-tof0
+    femtoReader->AddTrigId(370361); // tofmult3-vpd
+    femtoReader->AddTrigId(370341); // tofmult4
+
+    //
+    // Turn on only min-bias triggers
+    //
+    femtoReader->SetTrigger(370011, 1);
+    femtoReader->SetTrigger(370001, 1);
+
+    Long64_t mNEvents = femtoReader->GetNEvents();
+    mHbtManager->SetEventReader(femtoReader);
+
+    //
+    // Set basic event cut
+    //
+    gregsEventCut *mBasicEventCut = new gregsEventCut();
+    mBasicEventCut->SetCollisionType(0);     //0-pp, 1-AuAu
+    mBasicEventCut->SetEventMult(0, 45);
+    mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+    mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+    mBasicEventCut->SetVertRPos(mPrimVertR);
+
+    //
+    // Set basic kaon single-particle cut
+    //
+    gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+    mBasicTrackCut->SetCharge(mPosCharge);
+    mBasicTrackCut->SetMass(mChargedKaonMass);
+    mBasicTrackCut->SetType(mTrackType);
+    mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+    mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+    mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+    mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+    mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+    mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+    mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+    mBasicTrackCut->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+    mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+    mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+    mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+    mBasicTrackCut->SetTofMassSqr(mKaonTofMsqr[0],mKaonTofMsqr[1]);
+    mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+    mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+    mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+    mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+    mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+    mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+    mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+    //
+    // Set basic kaon pair cut
+    //
+    gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+    mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+    mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+    mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+    mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+    mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+    Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix, mCurrentPid;
+
+    //
+    // The multiplicity analysis matrix is CHARGE_BINS x MULT_BINS, where
+    // the first element corresponds to the particle charge: (pos, neg),
+    // and the second one corresponds to the RefMult intervals as listed below:
+    // Mult: 0-45, 0-2, 3-4, 5-6, 7-8, 9-10, 11-45
+    //
+    // StHbtAnalysis, cut, CF and monitor declarations
+    //
+    StHbtVertexMultAnalysis *mKaonMultAnalysis[CHARGE_BINS][SPH_BINS];
+
+    //
+    // Cuts
+    //
+    gregsEventCut       *mKaonEventCut[CHARGE_BINS][SPH_BINS];
+    gregsSimpleTrackCut *mKaonTrackCut[CHARGE_BINS][SPH_BINS];
+    gregsTrackPairCut   *mKaonPairCut[CHARGE_BINS][SPH_BINS];
+
+    //
+    // CorrFctns
+    //
+    QinvCorrFctnKt                *mKaonQinvMixKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnOpposHemisphereKt *mKaonQinvOppKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnRotationKt        *mKaonQinvRotKtCF[CHARGE_BINS][SPH_BINS];
+
+#ifdef MONITORS
+    //
+    // Cut monitors
+    //
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorFail[CHARGE_BINS][SPH_BINS];
+#endif
+
+    //
+    // Analysis creation
+    //
+    float sphLo, sphHi;
+    mMultLow = 0;
+    mMultHi = 50;
+    mMultBins = 10;
+    mEvents2Mix = 10;
+    for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        //
+        // sphLo < sphericity <= sphHi
+        //
+        switch (iSph)
+        {
+            case 0:  sphLo = -0.1;     sphHi = 0.1;  break;
+            default: sphLo = iSph*0.1; sphHi = (iSph + 1)*0.1;
+        }
+
+        for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            switch(iCharge)
+            {
+                case 0:   mLocalCharge = mPosCharge;    break;
+                default:  mLocalCharge = mNegCharge;    break;
+            };
+
+            //
+            // Create event cut
+            //
+            mKaonEventCut[iCharge][iSph] = new gregsEventCut(*mBasicEventCut);
+            mKaonEventCut[iCharge][iSph]->SetEventMult(mMultLow, mMultHi);
+            mKaonEventCut[iCharge][iSph]->SetSphericityCut(sphLo, sphHi);
+#ifdef MONITORS
+            mFxtKaonEventCutMonitorPass[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonPass", Form("_%d_%d", iCharge, iSph));
+            mFxtKaonEventCutMonitorFail[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonFail", Form("_%d_%d", iCharge, iSph));
+            mKaonEventCut[iCharge][iSph]->AddCutMonitor(mFxtKaonEventCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonEventCutMonitorFail[iCharge][iSph]);
+#endif
+
+            //
+            //  Like-sign kaon analysis  
+            //
+            // create analysis
+            mKaonMultAnalysis[iCharge][iSph] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+                                                                           mPrimVertZ[0],
+                                                                           mPrimVertZ[1],
+                                                                           mMultBins,
+                                                                           mMultLow,
+                                                                           mMultHi);
+
+            mKaonMultAnalysis[iCharge][iSph]->SetEventCut(mKaonEventCut[iCharge][iSph]);
+
+            //
+            // Create single-track cuts and add them to the analysis
+            //
+            mKaonTrackCut[iCharge][iSph] = new gregsSimpleTrackCut(*mBasicTrackCut);
+            mKaonTrackCut[iCharge][iSph]->SetCharge(mLocalCharge);
+#ifdef MONITORS
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonPass_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonFail_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mKaonTrackCut[iCharge][iSph]->AddCutMonitor(mFxtKaonTrackCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonTrackCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetFirstParticleCut(mKaonTrackCut[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->SetSecondParticleCut(mKaonTrackCut[iCharge][iSph]);
+
+            //
+            // Create pair cut
+            //
+            mKaonPairCut[iCharge][iSph] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+            mFxtKaonPairCutMonitorPass[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonPass_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mFxtKaonPairCutMonitorFail[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonFail_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mKaonPairCut[iCharge][iSph]->AddCutMonitor(mFxtKaonPairCutMonitorPass[iCharge][iSph],
+                                                       mFxtKaonPairCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetPairCut(mKaonPairCut[iCharge][iSph]);
+
+            //
+            // Correlation function initialization
+            //
+            mKaonQinvMixKtCF[iCharge][iSph] =
+                new QinvCorrFctnKt(Form("hKaonQinvMixKt_%d_%d", iCharge, iSph),
+                                   mQinvBins, mQinvVal[0], mQinvVal[1],
+                                   mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+            mKaonQinvOppKtCF[iCharge][iSph] =
+                new QinvCorrFctnOpposHemisphereKt(Form("hKaonQinvOppKt_%d_%d", iCharge, iSph),
+                                                  mQinvBins, mQinvVal[0], mQinvVal[1],
+                                                  mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                                  mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+            mKaonQinvRotKtCF[iCharge][iSph] =
+                new QinvCorrFctnRotationKt(Form("hKaonQinvRotKt_%d_%d", iCharge, iSph),
+                                           mQinvBins, mQinvVal[0], mQinvVal[1],
+                                           mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                           mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvMixKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvOppKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvRotKtCF[iCharge][iSph]);
+
+            mKaonMultAnalysis[iCharge][iSph]->SetNumEventsToMix(mEvents2Mix); // set number of events to mix
+            mHbtManager->AddAnalysis(mKaonMultAnalysis[iCharge][iSph]); // add analysis to the manager
+        }
+    }
+
+    mChain->Init(); // StChain initialization
+
+    //
+    // Start processing
+    //
+    Int_t iRet = 0;
+    Int_t mNEventsProcessed = 0;
+    Float_t mPercentCounter = 0.05;
+    Float_t mProgress = 0.;
+
+    unsigned int mSubCounter = 1;              // time
+    time_t mStartTime, mStopTime, mDiffTime;   // time
+    float mFrac;                               // time
+    mStartTime = time(0);                      // time
+
+    // mNEvents = 15;
+    std::cout << "Total Events: " << mNEvents << std::endl;
+    for(Int_t iEvent=0; iEvent < mNEvents; iEvent++)
+    {
+        mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+        mNEventsProcessed++;
+
+        if(mProgress >= mPercentCounter)
+        {
+            mPercentCounter += 0.05;
+            mStopTime = time(0);
+            mDiffTime = difftime(mStopTime, mStartTime);
+            mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+            mSubCounter++;
+            std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+                << std::endl;
+        }
+
+        mChain->Clear();
+        iRet = mChain->Make(iEvent);
+        if(iRet!=0)
+        {
+            std::cout << "[ERROR] - iRet: " << iRet << std::endl;
+            break;
+        }
+    }
+
+    mChain->Finish();
+
+    TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+
+    //
+    // Write hbt histograms to the file
+    //
+    std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+    for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexZ()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->RefMult()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexZ()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->RefMult()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Nhits()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->P()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Pt()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsDedx()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Rapidity()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsInvBeta()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Nhits()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->P()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Pt()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsDedx()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Rapidity()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsInvBeta()->Write();
+
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->Kt()->Write();
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->AveSepVsQinv()->Write();
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->Kt()->Write();
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->AveSepVsQinv()->Write();
+        }
+    }
+
+#endif
+
+    for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        for (Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            for (Int_t iKt = 0; iKt < mPairKtBins; iKt++)
+            {
+                mKaonQinvMixKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvMixKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mKaonQinvOppKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvOppKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mKaonQinvRotKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvRotKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+            }
+        }
+    }
+
+    std::cout << "\t[DONE]" << std::endl;
+
+    //
+    // Write the output file and close it
+    //
+    std::cout << "Saving the output file...";
+    oFile->Close();
+    std::cout << "\t[DONE]" << std::endl;
+
+    delete mChain;
+
+    std::cout << "********************************" << std::endl
+        << "*       expKaonHBT_pp200_1D    *" << std::endl
+        << "*              Finish          *" << std::endl
+        << "********************************" << std::endl << std::endl;
+}

+ 519 - 0
macros/HBT/pp/nomult/expKaonHBT_pp510_y2011_1D.C

@@ -0,0 +1,519 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+//
+// To define cut monitors
+// uncomments line below
+//
+#define MONITORS
+
+//
+// Constants for analyses
+//
+#define CHARGE_BINS 2
+#define SPH_BINS 10
+
+//
+//  Event information
+//
+const int mPrimVertZbins = 6;                // number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};    // redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                // distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};  // VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const bool mUseZdcCorrection = true;         // valid for AuAu collisions. True for 200GeV only!
+
+//
+// Track information
+//
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1;                     // primary track
+const int mPythiaKaonCode = 321;              // pythia particle code
+const int mPythiaPionCode = 211;              // pythia particle code
+
+const int mTrackNHits[2] = {10, 50};          // number of track hits
+const double mTrackNHitsFit[2] = {5, 50};    // number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};         // DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.;
+const double mTrackMom[2] = {0.12, 1.55};        // general momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.05, 1.55};  // general transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -0.5, 0.5};
+
+const short mDetectorSelection = 4;  // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  // track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2] = {-1e9, 1e9};
+const double mnSigmaKaon[2] = {-2., 2.};
+const double mnSigmaProton[2] = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.45};   // track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; // based on p+p data
+const double mKaonTofMsqr[2] = {0.2, 0.32};    // based on p+p data
+
+const double mTrackTnTMom[2] = {0.15, 1.55};
+const double mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-1e9, 1e9};
+const double mTnTNSigmaKaon[2] = {-2., 2.};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.5; // TPC&&TOF (P>Pthr)||TPC(P<Pthr)
+
+//
+// Pair information
+//
+const double mFmrCut[2] = {-1.1, 0.1};         // not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6};  // splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {-1e9., 1e9};        // pair transverse momentum cut
+const double mPairQinv[2] = {-1e9., 1e9};      // pair relative momentum
+const double mAveSep[2] = {5., 5000.};         // average separation of tracks within TPC
+
+//
+// Histograms values
+//
+const int mPairKtBins = 28;                // bin width 
+const double mPairKtVal[2] = {0.1, 1.5};   // 50 MeV/c 
+const int mQinvBins = 150;                 // 20 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 3.0};      // redefine for each analysis
+
+const Char_t *defaultInFile = "test.list";
+const Char_t *defaultOutFile = "oExpKaon_pp510_1D.root";
+
+
+
+void expKaonHBT_pp510_y2011_1D(const Char_t *inFileList = defaultInFile,
+                               const Char_t *oFileName  = defaultOutFile)
+{
+    std::cout << "*****************************" << std::endl
+              << "*    expKaonHBT_pp510_1D    *" << std::endl
+              << "*          Start            *" << std::endl
+              << "*****************************" << std::endl << std::endl;
+
+    //
+    // Load libraries
+    //
+    std::cout << "Loading libraries..." << std::endl;
+    gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+    loadSharedLibraries();
+    gSystem->Load("libMinuit");
+    gSystem->Load("StRefMultCorr");
+    gSystem->Load("StFlowMaker");   // should be placed before HBT
+    gSystem->Load("StHbtMaker");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("StPicoDstMakerRun12");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libStDb_Tables.so");
+    gSystem->Load("libgen_Tables.so");
+    gSystem->Load("libgeometry_Tables.so");
+    gSystem->Load("libsim_Tables.so");
+    gSystem->Load("libStarMagField.so");
+    gSystem->Load("libSt_db_Maker.so");
+    gSystem->Load("libSt_g2t.so");
+    gSystem->Load("libSt_geant_Maker.so");
+    gSystem->Load("StarGeneratorUtil.so");
+    gSystem->Load("StarGeneratorEvent.so");
+    gSystem->Load("StarGeneratorBase.so");
+    gSystem->Load("StFemtoDstMakerSph");
+
+    std::cout << "Libraries have been successfully loaded" << std::endl;
+
+    //
+    // Create chain
+    //
+    StChain* mChain = new StChain("StChain");
+    mChain->SetDebug(0);
+    StMuDebug::setLevel(0);
+
+    //
+    // Set StHbtMaker and StHbtManager
+    //
+    StHbtMaker *mHbtMaker = new StHbtMaker();
+    mHbtMaker->SetDebug(0);
+    StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+    //
+    // Set StHbtReader
+    //
+    StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+    femtoReader->AddTrigId(320000);   //Trigger selection for pp500 y2011
+    femtoReader->AddTrigId(320001);
+    femtoReader->AddTrigId(320011);
+    femtoReader->AddTrigId(320021);
+    femtoReader->AddTrigId(330021);
+
+    femtoReader->AddTrigId(320103);  // NOT VPDMB here and below
+    femtoReader->AddTrigId(320113);
+    femtoReader->AddTrigId(320123);
+    femtoReader->AddTrigId(330123);
+    femtoReader->AddTrigId(320500);
+    femtoReader->AddTrigId(320504);
+    femtoReader->AddTrigId(320514);
+    femtoReader->AddTrigId(320524);
+    femtoReader->AddTrigId(330524);
+    femtoReader->AddTrigId(320501);
+    femtoReader->AddTrigId(330501);
+    femtoReader->AddTrigId(320503);
+    femtoReader->AddTrigId(330503);
+
+    femtoReader->SetTrigger(320000, 1);   // set only VPDMB triggers
+    femtoReader->SetTrigger(320001, 1);
+    femtoReader->SetTrigger(320011, 1);
+    femtoReader->SetTrigger(320021, 1);
+    femtoReader->SetTrigger(330021, 1);
+
+    Long64_t mNEvents = femtoReader->GetNEvents();
+    mHbtManager->SetEventReader(femtoReader);
+
+    //
+    // Set basic event cut
+    //
+    gregsEventCut *mBasicEventCut = new gregsEventCut();
+    mBasicEventCut->SetCollisionType(0);     //0-pp, 1-AuAu
+    mBasicEventCut->SetEventMult(0, 45);
+    mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+    mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+    mBasicEventCut->SetVertRPos(mPrimVertR);
+
+    //
+    // Set basic kaon single-particle cut
+    //
+    gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+    mBasicTrackCut->SetCharge(mPosCharge);
+    mBasicTrackCut->SetMass(mChargedKaonMass);
+    mBasicTrackCut->SetType(mTrackType);
+    mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+    mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+    mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+    mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+    mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+    mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+    mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+    mBasicTrackCut->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+    mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+    mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+    mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+    mBasicTrackCut->SetTofMassSqr(mKaonTofMsqr[0],mKaonTofMsqr[1]);
+    mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+    mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+    mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+    mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+    mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+    mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+    mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+    //
+    // Set basic kaon pair cut
+    //
+    gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+    mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+    mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+    mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+    mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+    mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+    Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix, mCurrentPid;
+
+    //
+    // The multiplicity analysis matrix is CHARGE_BINS x MULT_BINS, where
+    // the first element corresponds to the particle charge: (pos, neg),
+    // and the second one corresponds to the RefMult intervals as listed below:
+    // Mult: 0-45, 0-2, 3-4, 5-6, 7-8, 9-10, 11-45
+    //
+    // StHbtAnalysis, cut, CF and monitor declarations
+    //
+    StHbtVertexMultAnalysis *mKaonMultAnalysis[CHARGE_BINS][SPH_BINS];
+
+    //
+    // Cuts
+    //
+    gregsEventCut       *mKaonEventCut[CHARGE_BINS][SPH_BINS];
+    gregsSimpleTrackCut *mKaonTrackCut[CHARGE_BINS][SPH_BINS];
+    gregsTrackPairCut   *mKaonPairCut[CHARGE_BINS][SPH_BINS];
+
+    //
+    // CorrFctns
+    //
+    QinvCorrFctnKt                *mKaonQinvMixKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnOpposHemisphereKt *mKaonQinvOppKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnRotationKt        *mKaonQinvRotKtCF[CHARGE_BINS][SPH_BINS];
+
+#ifdef MONITORS
+    //
+    // Cut monitors
+    //
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorFail[CHARGE_BINS][SPH_BINS];
+#endif
+
+    //
+    // Analysis creation
+    //
+    float sphLo, sphHi;
+    mMultLow = 0;
+    mMultHi = 50;
+    mMultBins = 10;
+    mEvents2Mix = 10;
+    for(Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        //
+        // sphLo < sphericity <= sphHi
+        //
+        switch (iSph)
+        {
+            case 0:  sphLo = -0.1;     sphHi = 0.1;  break;
+            default: sphLo = iSph*0.1; sphHi = (iSph + 1)*0.1;
+        }
+
+        for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            switch(iCharge)
+            {
+                case 0:   mLocalCharge = mPosCharge;    break;
+                default:  mLocalCharge = mNegCharge;    break;
+            };
+
+            //
+            // Create event cut
+            //
+            mKaonEventCut[iCharge][iSph] = new gregsEventCut(*mBasicEventCut);
+            mKaonEventCut[iCharge][iSph]->SetEventMult(mMultLow, mMultHi);
+            mKaonEventCut[iCharge][iSph]->SetSphericityCut(sphLo, sphHi);
+#ifdef MONITORS
+            mFxtKaonEventCutMonitorPass[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonPass", Form("_%d_%d", iCharge, iSph));
+            mFxtKaonEventCutMonitorFail[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonFail", Form("_%d_%d", iCharge, iSph));
+            mKaonEventCut[iCharge][iSph]->AddCutMonitor(mFxtKaonEventCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonEventCutMonitorFail[iCharge][iSph]);
+#endif
+            //
+            //  Like-sign kaon analysis
+            //
+            // create analysis
+            mKaonMultAnalysis[iCharge][iSph] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+                                                                           mPrimVertZ[0],
+                                                                           mPrimVertZ[1],
+                                                                           mMultBins,
+                                                                           mMultLow,
+                                                                           mMultHi);
+
+            mKaonMultAnalysis[iCharge][iSph]->SetEventCut(mKaonEventCut[iCharge][iSph]);
+
+            //
+            // Create single-track cuts and add them to the analysis
+            //
+            mKaonTrackCut[iCharge][iSph] = new gregsSimpleTrackCut(*mBasicTrackCut);
+            mKaonTrackCut[iCharge][iSph]->SetCharge(mLocalCharge);
+#ifdef MONITORS
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonPass_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonFail_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mKaonTrackCut[iCharge][iSph]->AddCutMonitor(mFxtKaonTrackCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonTrackCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetFirstParticleCut(mKaonTrackCut[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->SetSecondParticleCut(mKaonTrackCut[iCharge][iSph]);
+
+            //
+            // Create pair cut
+            //
+            mKaonPairCut[iCharge][iSph] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+            mFxtKaonPairCutMonitorPass[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonPass_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mFxtKaonPairCutMonitorFail[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonFail_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mKaonPairCut[iCharge][iSph]->AddCutMonitor(mFxtKaonPairCutMonitorPass[iCharge][iSph],
+                                                       mFxtKaonPairCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetPairCut(mKaonPairCut[iCharge][iSph]);
+
+            //
+            // Correlation function initialization
+            //
+            mKaonQinvMixKtCF[iCharge][iSph] =
+                new QinvCorrFctnKt(Form("hKaonQinvMixKt_%d_%d", iCharge, iSph),
+                                   mQinvBins, mQinvVal[0], mQinvVal[1],
+                                   mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+            mKaonQinvOppKtCF[iCharge][iSph] =
+                new QinvCorrFctnOpposHemisphereKt(Form("hKaonQinvOppKt_%d_%d", iCharge, iSph),
+                                                  mQinvBins, mQinvVal[0], mQinvVal[1],
+                                                  mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                                  mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+            mKaonQinvRotKtCF[iCharge][iSph] =
+                new QinvCorrFctnRotationKt(Form("hKaonQinvRotKt_%d_%d", iCharge, iSph),
+                                           mQinvBins, mQinvVal[0], mQinvVal[1],
+                                           mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                           mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvMixKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvOppKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvRotKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->SetNumEventsToMix(mEvents2Mix); // set number of events to mix
+            mHbtManager->AddAnalysis(mKaonMultAnalysis[iCharge][iSph]); // add analysis to the manager
+        }
+    }
+
+    mChain->Init(); // StChain initialization
+
+    //
+    // Start processing
+    //
+    Int_t iRet = 0;
+    Int_t mNEventsProcessed = 0;
+    Float_t mPercentCounter = 0.05;
+    Float_t mProgress = 0.;
+
+    unsigned int mSubCounter = 1;              // time
+    time_t mStartTime, mStopTime, mDiffTime;   // time
+    float mFrac;                               // time
+    mStartTime = time(0);                      // time
+
+    // mNEvents = 15;
+    std::cout << "Total Events: " << mNEvents << std::endl;
+    for(Int_t iEvent=0; iEvent < mNEvents; iEvent++)
+    {
+
+        mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+        mNEventsProcessed++;
+
+        if(mProgress >= mPercentCounter)
+        {
+            mPercentCounter += 0.05;
+            mStopTime = time(0);
+            mDiffTime = difftime(mStopTime, mStartTime);
+            mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+            mSubCounter++;
+            std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+                << std::endl;
+        }
+
+        mChain->Clear();
+        iRet = mChain->Make(iEvent);
+        if(iRet!=0)
+        {
+            std::cout << "[ERROR] - iRet: " << iRet << std::endl;
+            break;
+        }
+    }
+
+    mChain->Finish();
+
+    TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+
+    //
+    // Write hbt histograms to the file
+    //
+    std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+    for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+            for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+            {
+                mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexZ()->Write();
+                mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexYvsVertexX()->Write();
+                mFxtKaonEventCutMonitorPass[iCharge][iSph]->RefMult()->Write();
+                mFxtKaonEventCutMonitorPass[iCharge][iSph]->VpdVzDiff()->Write();
+                mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexZ()->Write();
+                mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexYvsVertexX()->Write();
+                mFxtKaonEventCutMonitorFail[iCharge][iSph]->RefMult()->Write();
+                mFxtKaonEventCutMonitorFail[iCharge][iSph]->VpdVzDiff()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->DCAGlobal()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Nhits()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->P()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Pt()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaPion()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaProton()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsDedx()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Rapidity()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PseudoRapidity()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsMassSqr()->Write();
+                mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsInvBeta()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->DCAGlobal()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Nhits()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->P()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Pt()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaPion()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaProton()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsDedx()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Rapidity()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PseudoRapidity()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsMassSqr()->Write();
+                mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsInvBeta()->Write();
+
+                mFxtKaonPairCutMonitorPass[iCharge][iSph]->Kt()->Write();
+                mFxtKaonPairCutMonitorPass[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+                mFxtKaonPairCutMonitorPass[iCharge][iSph]->FMRvsQinv()->Write();
+                mFxtKaonPairCutMonitorPass[iCharge][iSph]->SLvsQinv()->Write();
+                mFxtKaonPairCutMonitorPass[iCharge][iSph]->AveSepVsQinv()->Write();
+                mFxtKaonPairCutMonitorFail[iCharge][iSph]->Kt()->Write();
+                mFxtKaonPairCutMonitorFail[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+                mFxtKaonPairCutMonitorFail[iCharge][iSph]->FMRvsQinv()->Write();
+                mFxtKaonPairCutMonitorFail[iCharge][iSph]->SLvsQinv()->Write();
+                mFxtKaonPairCutMonitorFail[iCharge][iSph]->AveSepVsQinv()->Write();
+            }
+    }
+
+#endif
+    for (Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            for (Int_t iKt = 0; iKt < mPairKtBins; iKt++)
+            {
+                mKaonQinvMixKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvMixKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mKaonQinvOppKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvOppKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mKaonQinvRotKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mKaonQinvRotKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+            }
+        }
+    }
+
+    std::cout << "\t[DONE]" << std::endl;
+
+    //
+    // Write the output file and close it
+    //
+    std::cout << "Saving the output file...";
+    oFile->Close();
+    std::cout << "\t[DONE]" << std::endl;
+
+    delete mChain;
+
+    std::cout << "********************************" << std::endl
+        << "*       expKaonHBT_pp500_1D    *" << std::endl
+        << "*              Finish          *" << std::endl
+        << "********************************" << std::endl << std::endl;
+}

+ 519 - 0
macros/HBT/pp/nomult/expPionHBT_pp200_y2012_1D.C

@@ -0,0 +1,519 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+//
+// To define cut monitors
+// uncomments line below
+//
+#define MONITORS
+
+//
+// Constants for analyses
+//
+
+#define CHARGE_BINS 2
+#define SPH_BINS 10
+
+//
+//  Event information
+//
+const int mPrimVertZbins = 6;                // number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};    // redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                // distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};  // VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const bool mUseZdcCorrection = true;         // valid for AuAu collisions. True for 200GeV only!
+
+//
+// Track information
+//
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1;                     // primary track
+const int mPythiaKaonCode = 321;              // pythia particle code
+const int mPythiaPionCode = 211;              // pythia particle code
+
+const int mTrackNHits[2] = {10, 50};          // number of track hits
+const double mTrackNHitsFit[2] = {5, 50};    // number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};         // DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.;
+const double mTrackMom[2] = {0.12, 1.55};        // general momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.05, 1.55};  // general transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -0.5, 0.5};
+
+const short mDetectorSelection = 4;  // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  // track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2] = {-2., 2.};
+const double mnSigmaKaon[2] = {-1e9, 1e9};
+const double mnSigmaProton[2] = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.45};   // track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; // based on p+p data
+const double mKaonTofMsqr[2] = {0.2, 0.32};    // based on p+p data
+
+const double mTrackTnTMom[2] = {0.15, 1.55};
+const double mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-2., 2.};
+const double mTnTNSigmaKaon[2] = {-1e9, 1e9};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.5; // TPC&&TOF (P>Pthr)||TPC(P<Pthr)
+
+//
+// Pair information
+//
+const double mFmrCut[2] = {-1.1, 0.1};         // not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6};  // splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {-1e9., 1e9};        // pair transverse momentum cut
+const double mPairQinv[2] = {-1e9., 1e9};      // pair relative momentum
+const double mAveSep[2] = {5., 5000.};         // average separation of tracks within TPC
+
+//
+// Histograms values
+//
+const int mPairKtBins = 28;                // bin width 
+const double mPairKtVal[2] = {0.1, 1.5};   // 50 MeV/c 
+const int mQinvBins = 150;                 // 20 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 3.0};      // redefine for each analysis
+
+const Char_t *defaultInFile = "test.list";
+const Char_t *defaultOutFile = "oExpPion_pp200_1D.root";
+
+
+
+void expPionHBT_pp200_y2012_1D(const Char_t *inFileList = defaultInFile,
+                               const Char_t *oFileName  = defaultOutFile)
+{
+    std::cout << "*****************************" << std::endl
+              << "*    expPionHBT_pp200_1D    *" << std::endl
+              << "*          Start            *" << std::endl
+              << "*****************************" << std::endl << std::endl;
+
+    //
+    // Load libraries
+    //
+    std::cout << "Loading libraries..." << std::endl;
+    gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+    loadSharedLibraries();
+    gSystem->Load("libMinuit");
+    gSystem->Load("StRefMultCorr");
+    gSystem->Load("StFlowMaker");   // should be placed before HBT
+    gSystem->Load("StHbtMaker");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("StPicoDstMakerRun12");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libStDb_Tables.so");
+    gSystem->Load("libgen_Tables.so");
+    gSystem->Load("libgeometry_Tables.so");
+    gSystem->Load("libsim_Tables.so");
+    gSystem->Load("libStarMagField.so");
+    gSystem->Load("libSt_db_Maker.so");
+    gSystem->Load("libSt_g2t.so");
+    gSystem->Load("libSt_geant_Maker.so");
+    gSystem->Load("StarGeneratorUtil.so");
+    gSystem->Load("StarGeneratorEvent.so");
+    gSystem->Load("StarGeneratorBase.so");
+    gSystem->Load("StFemtoDstMakerSph");
+
+    std::cout << "Libraries have been successfully loaded" << std::endl;
+
+    //
+    // Create chain
+    //
+    StChain* mChain = new StChain("StChain");
+    mChain->SetDebug(0);
+    StMuDebug::setLevel(0);
+
+    //
+    // Set StHbtMaker and StHbtManager
+    //
+    StHbtMaker *mHbtMaker = new StHbtMaker();
+    mHbtMaker->SetDebug(0);
+    StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+    //
+    // Set StHbtReader
+    //
+    StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+    femtoReader->AddTrigId(370011); // vpdmb-nobsmd
+    femtoReader->AddTrigId(370001); // vpdmb
+    femtoReader->AddTrigId(370022); // bbcmb
+    femtoReader->AddTrigId(370501); // bht0-vpdmb
+    femtoReader->AddTrigId(370511); // bht1-vpdmb
+    femtoReader->AddTrigId(370531); // bht2
+    femtoReader->AddTrigId(370542); // bht0-bbcmb-tof0
+    femtoReader->AddTrigId(370546); // bht1-bbcmb-tof0
+    femtoReader->AddTrigId(370522); // bht2-bbcmb
+    femtoReader->AddTrigId(370601); // jp0
+    femtoReader->AddTrigId(370611); // jp1
+    femtoReader->AddTrigId(370621); // jp2
+    femtoReader->AddTrigId(370641); // ajp
+    femtoReader->AddTrigId(370301); // bbcmb-tof0
+    femtoReader->AddTrigId(370361); // tofmult3-vpd
+    femtoReader->AddTrigId(370341); // tofmult4
+
+    //
+    // Turn on only min-bias triggers
+    //
+    femtoReader->SetTrigger(370011, 1);
+    femtoReader->SetTrigger(370001, 1);
+
+    Long64_t mNEvents = femtoReader->GetNEvents();
+    mHbtManager->SetEventReader(femtoReader);
+
+    //
+    // Set basic event cut
+    //
+    gregsEventCut *mBasicEventCut = new gregsEventCut();
+    mBasicEventCut->SetCollisionType(0);     //0-pp, 1-AuAu
+    mBasicEventCut->SetEventMult(0, 45);
+    mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+    mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+    mBasicEventCut->SetVertRPos(mPrimVertR);
+
+    //
+    // Set basic kaon single-particle cut
+    //
+    gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+    mBasicTrackCut->SetCharge(mPosCharge);
+    mBasicTrackCut->SetMass(mChargedPionMass);
+    mBasicTrackCut->SetType(mTrackType);
+    mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+    mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+    mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+    mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+    mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+    mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+    mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+    mBasicTrackCut->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+    mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+    mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+    mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+    mBasicTrackCut->SetTofMassSqr(mPionTofMsqr[0],mPionTofMsqr[1]);
+    mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+    mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+    mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+    mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+    mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+    mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+    mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+    //
+    // Set basic kaon pair cut
+    //
+    gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+    mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+    mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+    mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+    mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+    mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+    Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix, mCurrentPid;
+
+    //
+    // The multiplicity analysis matrix is CHARGE_BINS x MULT_BINS, where
+    // the first element corresponds to the particle charge: (pos, neg),
+    // and the second one corresponds to the RefMult intervals as listed below:
+    // Mult: 0-45, 0-2, 3-4, 5-6, 7-8, 9-10, 11-45
+    //
+    // StHbtAnalysis, cut, CF and monitor declarations
+    //
+    StHbtVertexMultAnalysis *mPionMultAnalysis[CHARGE_BINS][SPH_BINS];
+
+    //
+    // Cuts
+    //
+    gregsEventCut       *mPionEventCut[CHARGE_BINS][SPH_BINS];
+    gregsSimpleTrackCut *mPionTrackCut[CHARGE_BINS][SPH_BINS];
+    gregsTrackPairCut   *mPionPairCut[CHARGE_BINS][SPH_BINS];
+
+    //
+    // CorrFctns
+    //
+    QinvCorrFctnKt                *mPionQinvMixKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnOpposHemisphereKt *mPionQinvOppKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnRotationKt        *mPionQinvRotKtCF[CHARGE_BINS][SPH_BINS];
+
+#ifdef MONITORS
+    //
+    // Cut monitors
+    //
+    fxtEventCutMonitor *mFxtPionEventCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtPionTrackCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtPionPairCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtEventCutMonitor *mFxtPionEventCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtPionTrackCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtPionPairCutMonitorFail[CHARGE_BINS][SPH_BINS];
+#endif
+
+    //
+    // Analysis creation
+    //
+    float sphLo, sphHi;
+    mMultLow = 0;
+    mMultHi = 50;
+    mMultBins = 10;
+    mEvents2Mix = 10;
+    for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        //
+        // sphLo < sphericity <= sphHi
+        //
+        switch (iSph)
+        {
+            case 0:  sphLo = -0.1;     sphHi = 0.1;  break;
+            default: sphLo = iSph*0.1; sphHi = (iSph + 1)*0.1;
+        }
+
+        for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            switch(iCharge)
+            {
+                case 0:   mLocalCharge = mPosCharge;    break;
+                default:  mLocalCharge = mNegCharge;    break;
+            };
+
+            //
+            // Create event cut
+            //
+            mPionEventCut[iCharge][iSph] = new gregsEventCut(*mBasicEventCut);
+            mPionEventCut[iCharge][iSph]->SetEventMult(mMultLow, mMultHi);
+            mPionEventCut[iCharge][iSph]->SetSphericityCut(sphLo, sphHi);
+#ifdef MONITORS
+            mFxtPionEventCutMonitorPass[iCharge][iSph] =
+                new fxtEventCutMonitor("eventPionPass", Form("_%d_%d", iCharge, iSph));
+            mFxtPionEventCutMonitorFail[iCharge][iSph] =
+                new fxtEventCutMonitor("eventPionFail", Form("_%d_%d", iCharge, iSph));
+            mPionEventCut[iCharge][iSph]->AddCutMonitor(mFxtPionEventCutMonitorPass[iCharge][iSph],
+                                                        mFxtPionEventCutMonitorFail[iCharge][iSph]);
+#endif
+
+            //
+            //  Like-sign kaon analysis  
+            //
+            // create analysis
+            mPionMultAnalysis[iCharge][iSph] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+                                                                           mPrimVertZ[0],
+                                                                           mPrimVertZ[1],
+                                                                           mMultBins,
+                                                                           mMultLow,
+                                                                           mMultHi);
+
+            mPionMultAnalysis[iCharge][iSph]->SetEventCut(mPionEventCut[iCharge][iSph]);
+
+            //
+            // Create single-track cuts and add them to the analysis
+            //
+            mPionTrackCut[iCharge][iSph] = new gregsSimpleTrackCut(*mBasicTrackCut);
+            mPionTrackCut[iCharge][iSph]->SetCharge(mLocalCharge);
+#ifdef MONITORS
+            mFxtPionTrackCutMonitorPass[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackPionPass_%d_%d_", iCharge, iSph), mChargedPionMass);
+            mFxtPionTrackCutMonitorFail[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackPionFail_%d_%d_", iCharge, iSph), mChargedPionMass);
+            mPionTrackCut[iCharge][iSph]->AddCutMonitor(mFxtPionTrackCutMonitorPass[iCharge][iSph],
+                                                        mFxtPionTrackCutMonitorFail[iCharge][iSph]);
+#endif
+            mPionMultAnalysis[iCharge][iSph]->SetFirstParticleCut(mPionTrackCut[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->SetSecondParticleCut(mPionTrackCut[iCharge][iSph]);
+
+            //
+            // Create pair cut
+            //
+            mPionPairCut[iCharge][iSph] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+            mFxtPionPairCutMonitorPass[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairPionPass_%d_%d_", iCharge, iSph));
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->SetParticleMass(mChargedPionMass);
+            mFxtPionPairCutMonitorFail[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairPionFail_%d_%d_", iCharge, iSph));
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->SetParticleMass(mChargedPionMass);
+            mPionPairCut[iCharge][iSph]->AddCutMonitor(mFxtPionPairCutMonitorPass[iCharge][iSph],
+                                                       mFxtPionPairCutMonitorFail[iCharge][iSph]);
+#endif
+            mPionMultAnalysis[iCharge][iSph]->SetPairCut(mPionPairCut[iCharge][iSph]);
+
+            //
+            // Correlation function initialization
+            //
+            mPionQinvMixKtCF[iCharge][iSph] =
+                new QinvCorrFctnKt(Form("hPionQinvMixKt_%d_%d", iCharge, iSph),
+                                   mQinvBins, mQinvVal[0], mQinvVal[1],
+                                   mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+            mPionQinvOppKtCF[iCharge][iSph] =
+                new QinvCorrFctnOpposHemisphereKt(Form("hPionQinvOppKt_%d_%d", iCharge, iSph),
+                                                  mQinvBins, mQinvVal[0], mQinvVal[1],
+                                                  mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                                  mChargedPionMass, mPionPairCut[iCharge][iSph]);
+            mPionQinvRotKtCF[iCharge][iSph] =
+                new QinvCorrFctnRotationKt(Form("hPionQinvRotKt_%d_%d", iCharge, iSph),
+                                           mQinvBins, mQinvVal[0], mQinvVal[1],
+                                           mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                           mChargedPionMass, mPionPairCut[iCharge][iSph]);
+
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvMixKtCF[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvOppKtCF[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvRotKtCF[iCharge][iSph]);
+
+            mPionMultAnalysis[iCharge][iSph]->SetNumEventsToMix(mEvents2Mix); // set number of events to mix
+            mHbtManager->AddAnalysis(mPionMultAnalysis[iCharge][iSph]); // add analysis to the manager
+        }
+    }
+
+    mChain->Init(); // StChain initialization
+
+    //
+    // Start processing
+    //
+    Int_t iRet = 0;
+    Int_t mNEventsProcessed = 0;
+    Float_t mPercentCounter = 0.05;
+    Float_t mProgress = 0.;
+
+    unsigned int mSubCounter = 1;              // time
+    time_t mStartTime, mStopTime, mDiffTime;   // time
+    float mFrac;                               // time
+    mStartTime = time(0);                      // time
+
+    // mNEvents = 15;
+    std::cout << "Total Events: " << mNEvents << std::endl;
+    for(Int_t iEvent=0; iEvent < mNEvents; iEvent++)
+    {
+        mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+        mNEventsProcessed++;
+
+        if(mProgress >= mPercentCounter)
+        {
+            mPercentCounter += 0.05;
+            mStopTime = time(0);
+            mDiffTime = difftime(mStopTime, mStartTime);
+            mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+            mSubCounter++;
+            std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+                << std::endl;
+        }
+
+        mChain->Clear();
+        iRet = mChain->Make(iEvent);
+        if(iRet!=0)
+        {
+            std::cout << "[ERROR] - iRet: " << iRet << std::endl;
+            break;
+        }
+    }
+
+    mChain->Finish();
+
+    TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+
+    //
+    // Write hbt histograms to the file
+    //
+    std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+    for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VertexZ()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->RefMult()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VertexZ()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->RefMult()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Nhits()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->P()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Pt()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsDedx()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Rapidity()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsInvBeta()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Nhits()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->P()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Pt()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsDedx()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Rapidity()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsInvBeta()->Write();
+
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->Kt()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->AveSepVsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->Kt()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->AveSepVsQinv()->Write();
+        }
+    }
+
+#endif
+
+    for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        for (Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            for (Int_t iKt = 0; iKt < mPairKtBins; iKt++)
+            {
+                mPionQinvMixKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvMixKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mPionQinvOppKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvOppKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mPionQinvRotKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvRotKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+            }
+        }
+    }
+
+    std::cout << "\t[DONE]" << std::endl;
+
+    //
+    // Write the output file and close it
+    //
+    std::cout << "Saving the output file...";
+    oFile->Close();
+    std::cout << "\t[DONE]" << std::endl;
+
+    delete mChain;
+
+    std::cout << "********************************" << std::endl
+        << "*       expPionHBT_pp200_1D    *" << std::endl
+        << "*              Finish          *" << std::endl
+        << "********************************" << std::endl << std::endl;
+}

+ 521 - 0
macros/HBT/pp/nomult/expPionHBT_pp510_y2011_1D.C

@@ -0,0 +1,521 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+//
+// To define cut monitors
+// uncomments line below
+//
+#define MONITORS
+
+//
+// Constants for analyses
+//
+
+#define CHARGE_BINS 2
+#define SPH_BINS 10
+
+//
+//  Event information
+//
+const int mPrimVertZbins = 6;                // number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};    // redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                // distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};  // VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const bool mUseZdcCorrection = true;         // valid for AuAu collisions. True for 200GeV only!
+
+//
+// Track information
+//
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1;                     // primary track
+const int mPythiaKaonCode = 321;              // pythia particle code
+const int mPythiaPionCode = 211;              // pythia particle code
+
+const int mTrackNHits[2] = {10, 50};          // number of track hits
+const double mTrackNHitsFit[2] = {5, 50};    // number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};         // DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.;
+const double mTrackMom[2] = {0.12, 1.55};        // general momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.05, 1.55};  // general transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -0.5, 0.5};
+
+const short mDetectorSelection = 4;  // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  // track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2] = {-2., 2.};
+const double mnSigmaKaon[2] = {-1e9, 1e9};
+const double mnSigmaProton[2] = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.45};   // track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; // based on p+p data
+const double mKaonTofMsqr[2] = {0.2, 0.32};    // based on p+p data
+
+const double mTrackTnTMom[2] = {0.15, 1.55};
+const double mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-2., 2.};
+const double mTnTNSigmaKaon[2] = {-1e9, 1e9};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.5; // TPC&&TOF (P>Pthr)||TPC(P<Pthr)
+
+//
+// Pair information
+//
+const double mFmrCut[2] = {-1.1, 0.1};         // not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6};  // splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {-1e9., 1e9};        // pair transverse momentum cut
+const double mPairQinv[2] = {-1e9., 1e9};      // pair relative momentum
+const double mAveSep[2] = {5., 5000.};         // average separation of tracks within TPC
+
+//
+// Histograms values
+//
+const int mPairKtBins = 28;                // bin width 
+const double mPairKtVal[2] = {0.1, 1.5};   // 50 MeV/c 
+const int mQinvBins = 150;                 // 20 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 3.0};      // redefine for each analysis
+
+const Char_t *defaultInFile = "test.list";
+const Char_t *defaultOutFile = "oExpPion_pp510_1D.root";
+
+
+
+void expPionHBT_pp510_y2011_1D(const Char_t *inFileList = defaultInFile,
+                               const Char_t *oFileName  = defaultOutFile)
+{
+    std::cout << "*****************************" << std::endl
+              << "*    expPionHBT_pp510_1D    *" << std::endl
+              << "*          Start            *" << std::endl
+              << "*****************************" << std::endl << std::endl;
+
+    //
+    // Load libraries
+    //
+    std::cout << "Loading libraries..." << std::endl;
+    gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+    loadSharedLibraries();
+    gSystem->Load("libMinuit");
+    gSystem->Load("StRefMultCorr");
+    gSystem->Load("StFlowMaker");   // should be placed before HBT
+    gSystem->Load("StHbtMaker");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("StPicoDstMakerRun12");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libStDb_Tables.so");
+    gSystem->Load("libgen_Tables.so");
+    gSystem->Load("libgeometry_Tables.so");
+    gSystem->Load("libsim_Tables.so");
+    gSystem->Load("libStarMagField.so");
+    gSystem->Load("libSt_db_Maker.so");
+    gSystem->Load("libSt_g2t.so");
+    gSystem->Load("libSt_geant_Maker.so");
+    gSystem->Load("StarGeneratorUtil.so");
+    gSystem->Load("StarGeneratorEvent.so");
+    gSystem->Load("StarGeneratorBase.so");
+    gSystem->Load("StFemtoDstMakerSph");
+
+    std::cout << "Libraries have been successfully loaded" << std::endl;
+
+    //
+    // Create chain
+    //
+    StChain* mChain = new StChain("StChain");
+    mChain->SetDebug(0);
+    StMuDebug::setLevel(0);
+
+    //
+    // Set StHbtMaker and StHbtManager
+    //
+    StHbtMaker *mHbtMaker = new StHbtMaker();
+    mHbtMaker->SetDebug(0);
+    StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+    //
+    // Set StHbtReader
+    //
+    StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+    femtoReader->AddTrigId(320000);   //Trigger selection for pp500 y2011
+    femtoReader->AddTrigId(320001);
+    femtoReader->AddTrigId(320011);
+    femtoReader->AddTrigId(320021);
+    femtoReader->AddTrigId(330021);
+
+    femtoReader->AddTrigId(320103);  // NOT VPDMB here and below
+    femtoReader->AddTrigId(320113);
+    femtoReader->AddTrigId(320123);
+    femtoReader->AddTrigId(330123);
+    femtoReader->AddTrigId(320500);
+    femtoReader->AddTrigId(320504);
+    femtoReader->AddTrigId(320514);
+    femtoReader->AddTrigId(320524);
+    femtoReader->AddTrigId(330524);
+    femtoReader->AddTrigId(320501);
+    femtoReader->AddTrigId(330501);
+    femtoReader->AddTrigId(320503);
+    femtoReader->AddTrigId(330503);
+
+    femtoReader->SetTrigger(320000, 1);   // set only VPDMB triggers
+    femtoReader->SetTrigger(320001, 1);
+    femtoReader->SetTrigger(320011, 1);
+    femtoReader->SetTrigger(320021, 1);
+    femtoReader->SetTrigger(330021, 1);
+
+    Long64_t mNEvents = femtoReader->GetNEvents();
+    mHbtManager->SetEventReader(femtoReader);
+
+    //
+    // Set basic event cut
+    //
+    gregsEventCut *mBasicEventCut = new gregsEventCut();
+    mBasicEventCut->SetCollisionType(0);     //0-pp, 1-AuAu
+    mBasicEventCut->SetEventMult(0, 45);
+    mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+    mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+    mBasicEventCut->SetVertRPos(mPrimVertR);
+
+    //
+    // Set basic kaon single-particle cut
+    //
+    gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+    mBasicTrackCut->SetCharge(mPosCharge);
+    mBasicTrackCut->SetMass(mChargedPionMass);
+    mBasicTrackCut->SetType(mTrackType);
+    mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+    mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+    mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+    mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+    mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+    mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+    mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+    mBasicTrackCut->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+    mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+    mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+    mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+    mBasicTrackCut->SetTofMassSqr(mPionTofMsqr[0],mPionTofMsqr[1]);
+    mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+    mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+    mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+    mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+    mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+    mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+    mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+    //
+    // Set basic kaon pair cut
+    //
+    gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+    mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+    mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+    mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+    mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+    mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+    Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix, mCurrentPid;
+
+    //
+    // The multiplicity analysis matrix is CHARGE_BINS x MULT_BINS, where
+    // the first element corresponds to the particle charge: (pos, neg),
+    // and the second one corresponds to the RefMult intervals as listed below:
+    // Mult: 0-45, 0-2, 3-4, 5-6, 7-8, 9-10, 11-45
+    //
+    // StHbtAnalysis, cut, CF and monitor declarations
+    //
+    StHbtVertexMultAnalysis *mPionMultAnalysis[CHARGE_BINS][SPH_BINS];
+
+    //
+    // Cuts
+    //
+    gregsEventCut       *mPionEventCut[CHARGE_BINS][SPH_BINS];
+    gregsSimpleTrackCut *mPionTrackCut[CHARGE_BINS][SPH_BINS];
+    gregsTrackPairCut   *mPionPairCut[CHARGE_BINS][SPH_BINS];
+
+    //
+    // CorrFctns
+    //
+    QinvCorrFctnKt                *mPionQinvMixKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnOpposHemisphereKt *mPionQinvOppKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnRotationKt        *mPionQinvRotKtCF[CHARGE_BINS][SPH_BINS];
+
+#ifdef MONITORS
+    //
+    // Cut monitors
+    //
+    fxtEventCutMonitor *mFxtPionEventCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtPionTrackCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtPionPairCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtEventCutMonitor *mFxtPionEventCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtPionTrackCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtPionPairCutMonitorFail[CHARGE_BINS][SPH_BINS];
+#endif
+
+    //
+    // Analysis creation
+    //
+    float sphLo, sphHi;
+    mMultLow = 0;
+    mMultHi = 50;
+    mMultBins = 10;
+    mEvents2Mix = 10;
+    for(Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        //
+        // sphLo < sphericity <= sphHi
+        //
+        switch (iSph)
+        {
+            case 0:  sphLo = -0.1;     sphHi = 0.1;  break;
+            default: sphLo = iSph*0.1; sphHi = (iSph + 1)*0.1;
+        }
+
+        for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            switch(iCharge)
+            {
+                case 0:   mLocalCharge = mPosCharge;    break;
+                default:  mLocalCharge = mNegCharge;    break;
+            };
+
+            //
+            // Create event cut
+            //
+            mPionEventCut[iCharge][iSph] = new gregsEventCut(*mBasicEventCut);
+            mPionEventCut[iCharge][iSph]->SetEventMult(mMultLow, mMultHi);
+            mPionEventCut[iCharge][iSph]->SetSphericityCut(sphLo, sphHi);
+#ifdef MONITORS
+            mFxtPionEventCutMonitorPass[iCharge][iSph] =
+                new fxtEventCutMonitor("eventPionPass", Form("_%d_%d", iCharge, iSph));
+            mFxtPionEventCutMonitorFail[iCharge][iSph] =
+                new fxtEventCutMonitor("eventPionFail", Form("_%d_%d", iCharge, iSph));
+            mPionEventCut[iCharge][iSph]->AddCutMonitor(mFxtPionEventCutMonitorPass[iCharge][iSph],
+                                                        mFxtPionEventCutMonitorFail[iCharge][iSph]);
+#endif
+            //
+            //  Like-sign kaon analysis
+            //
+            // create analysis
+            mPionMultAnalysis[iCharge][iSph] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+                                                                           mPrimVertZ[0],
+                                                                           mPrimVertZ[1],
+                                                                           mMultBins,
+                                                                           mMultLow,
+                                                                           mMultHi);
+
+            mPionMultAnalysis[iCharge][iSph]->SetEventCut(mPionEventCut[iCharge][iSph]);
+
+            //
+            // Create single-track cuts and add them to the analysis
+            //
+            mPionTrackCut[iCharge][iSph] = new gregsSimpleTrackCut(*mBasicTrackCut);
+            mPionTrackCut[iCharge][iSph]->SetCharge(mLocalCharge);
+#ifdef MONITORS
+            mFxtPionTrackCutMonitorPass[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackPionPass_%d_%d_", iCharge, iSph), mChargedPionMass);
+            mFxtPionTrackCutMonitorFail[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackPionFail_%d_%d_", iCharge, iSph), mChargedPionMass);
+            mPionTrackCut[iCharge][iSph]->AddCutMonitor(mFxtPionTrackCutMonitorPass[iCharge][iSph],
+                                                        mFxtPionTrackCutMonitorFail[iCharge][iSph]);
+#endif
+            mPionMultAnalysis[iCharge][iSph]->SetFirstParticleCut(mPionTrackCut[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->SetSecondParticleCut(mPionTrackCut[iCharge][iSph]);
+
+            //
+            // Create pair cut
+            //
+            mPionPairCut[iCharge][iSph] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+            mFxtPionPairCutMonitorPass[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairPionPass_%d_%d_", iCharge, iSph));
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->SetParticleMass(mChargedPionMass);
+            mFxtPionPairCutMonitorFail[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairPionFail_%d_%d_", iCharge, iSph));
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->SetParticleMass(mChargedPionMass);
+            mPionPairCut[iCharge][iSph]->AddCutMonitor(mFxtPionPairCutMonitorPass[iCharge][iSph],
+                                                       mFxtPionPairCutMonitorFail[iCharge][iSph]);
+#endif
+            mPionMultAnalysis[iCharge][iSph]->SetPairCut(mPionPairCut[iCharge][iSph]);
+
+            //
+            // Correlation function initialization
+            //
+            mPionQinvMixKtCF[iCharge][iSph] =
+                new QinvCorrFctnKt(Form("hPionQinvMixKt_%d_%d", iCharge, iSph),
+                                   mQinvBins, mQinvVal[0], mQinvVal[1],
+                                   mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+            mPionQinvOppKtCF[iCharge][iSph] =
+                new QinvCorrFctnOpposHemisphereKt(Form("hPionQinvOppKt_%d_%d", iCharge, iSph),
+                                                  mQinvBins, mQinvVal[0], mQinvVal[1],
+                                                  mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                                  mChargedPionMass, mPionPairCut[iCharge][iSph]);
+            mPionQinvRotKtCF[iCharge][iSph] =
+                new QinvCorrFctnRotationKt(Form("hPionQinvRotKt_%d_%d", iCharge, iSph),
+                                           mQinvBins, mQinvVal[0], mQinvVal[1],
+                                           mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                           mChargedPionMass, mPionPairCut[iCharge][iSph]);
+
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvMixKtCF[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvOppKtCF[iCharge][iSph]);
+            mPionMultAnalysis[iCharge][iSph]->AddCorrFctn(mPionQinvRotKtCF[iCharge][iSph]);
+
+            mPionMultAnalysis[iCharge][iSph]->SetNumEventsToMix(mEvents2Mix); // set number of events to mix
+            mHbtManager->AddAnalysis(mPionMultAnalysis[iCharge][iSph]); // add analysis to the manager
+        }
+    }
+
+    mChain->Init(); // StChain initialization
+
+    //
+    // Start processing
+    //
+    Int_t iRet = 0;
+    Int_t mNEventsProcessed = 0;
+    Float_t mPercentCounter = 0.05;
+    Float_t mProgress = 0.;
+
+    unsigned int mSubCounter = 1;              // time
+    time_t mStartTime, mStopTime, mDiffTime;   // time
+    float mFrac;                               // time
+    mStartTime = time(0);                      // time
+
+    // mNEvents = 15;
+    std::cout << "Total Events: " << mNEvents << std::endl;
+    for(Int_t iEvent=0; iEvent < mNEvents; iEvent++)
+    {
+
+        mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+        mNEventsProcessed++;
+
+        if(mProgress >= mPercentCounter)
+        {
+            mPercentCounter += 0.05;
+            mStopTime = time(0);
+            mDiffTime = difftime(mStopTime, mStartTime);
+            mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+            mSubCounter++;
+            std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+                << std::endl;
+        }
+
+        mChain->Clear();
+        iRet = mChain->Make(iEvent);
+        if(iRet!=0)
+        {
+            std::cout << "[ERROR] - iRet: " << iRet << std::endl;
+            break;
+        }
+    }
+
+    mChain->Finish();
+
+    TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+
+    //
+    // Write hbt histograms to the file
+    //
+    std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+    for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VertexZ()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->RefMult()->Write();
+            mFxtPionEventCutMonitorPass[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VertexZ()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->RefMult()->Write();
+            mFxtPionEventCutMonitorFail[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Nhits()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->P()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Pt()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsDedx()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->Rapidity()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtPionTrackCutMonitorPass[iCharge][iSph]->PvsInvBeta()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Nhits()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->P()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Pt()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsDedx()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->Rapidity()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtPionTrackCutMonitorFail[iCharge][iSph]->PvsInvBeta()->Write();
+
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->Kt()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtPionPairCutMonitorPass[iCharge][iSph]->AveSepVsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->Kt()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->Pt1Pt2Qinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->FMRvsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->SLvsQinv()->Write();
+            mFxtPionPairCutMonitorFail[iCharge][iSph]->AveSepVsQinv()->Write();
+        }
+    }
+
+#endif
+    for (Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            for (Int_t iKt = 0; iKt < mPairKtBins; iKt++)
+            {
+                mPionQinvMixKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvMixKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mPionQinvOppKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvOppKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+                mPionQinvRotKtCF[iCharge][iSph]->Numerator(iKt)->Write();
+                mPionQinvRotKtCF[iCharge][iSph]->Denominator(iKt)->Write();
+            }
+        }
+    }
+
+    std::cout << "\t[DONE]" << std::endl;
+
+    //
+    // Write the output file and close it
+    //
+    std::cout << "Saving the output file...";
+    oFile->Close();
+    std::cout << "\t[DONE]" << std::endl;
+
+    delete mChain;
+
+    std::cout << "********************************" << std::endl
+        << "*       expPionHBT_pp500_1D    *" << std::endl
+        << "*              Finish          *" << std::endl
+        << "********************************" << std::endl << std::endl;
+}

+ 497 - 0
macros/HBT/pp/nomult/simKaonHBT_pp200n510_1D.C

@@ -0,0 +1,497 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+//
+// To define cut monitors
+// uncomments line below
+//
+#define MONITORS
+
+//
+// Constants for analyses
+//
+
+#define CHARGE_BINS 2
+#define SPH_BINS 10
+
+//
+//  Event information
+//
+const int mPrimVertZbins = 6;                // number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};    // redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                // distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};  // VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const bool mUseZdcCorrection = true;         // valid for AuAu collisions. True for 200GeV only!
+
+//
+// Track information
+//
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1;                     // primary track
+const int mPythiaKaonCode = 321;              // pythia particle code
+const int mPythiaPionCode = 211;              // pythia particle code
+
+const int mTrackNHits[2] = {10, 50};          // number of track hits
+const double mTrackNHitsFit[2] = {5, 50};    // number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};         // DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.;
+const double mTrackMom[2] = {0.12, 1.55};        // general momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.05, 1.55};  // general transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -0.5, 0.5};
+
+const short mDetectorSelection = 4;  // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  // track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2] = {-1e9, 1e9};
+const double mnSigmaKaon[2] = {-2., 2.};
+const double mnSigmaProton[2] = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.45};   // track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; // based on p+p data
+const double mKaonTofMsqr[2] = {0.2, 0.32};    // based on p+p data
+
+const double mTrackTnTMom[2] = {0.15, 1.55};
+const double mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-1e9, 1e9};
+const double mTnTNSigmaKaon[2] = {-2., 2.};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.5; // TPC&&TOF (P>Pthr)||TPC(P<Pthr)
+
+//
+// Pair information
+//
+const double mFmrCut[2] = {-1.1, 0.1};         // not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6};  // splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {-1e9., 1e9};        // pair transverse momentum cut
+const double mPairQinv[2] = {-1e9., 1e9};      // pair relative momentum
+const double mAveSep[2] = {5., 5000.};         // average separation of tracks within TPC
+
+//
+// Histograms values
+//
+const int mPairKtBins = 28;                // bin width 
+const double mPairKtVal[2] = {0.1, 1.5};   // 50 MeV/c 
+const int mQinvBins = 150;                 // 20 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 3.0};      // redefine for each analysis
+
+const Char_t *defaultInFile = "test.list";
+const Char_t *defaultOutFile = "oSimKaon_pp200_1D.root";
+
+
+
+void simKaonHBT_pp200n510_1D(const Char_t *inFileList = defaultInFile,
+                             const Char_t *oFileName  = defaultOutFile)
+{
+    std::cout << "*************************************" << std::endl
+              << "*    simKaonHBT_pp200_and_510_1D    *" << std::endl
+              << "*                Start              *" << std::endl
+              << "*************************************" << std::endl << std::endl;
+
+    //
+    // Load libraries
+    //
+    std::cout << "Loading libraries..." << std::endl;
+    gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+    loadSharedLibraries();
+    gSystem->Load("libMinuit");
+    gSystem->Load("StRefMultCorr");
+    gSystem->Load("StFlowMaker");   // should be placed before HBT
+    gSystem->Load("StHbtMaker");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("StPicoDstMakerRun12");
+    gSystem->Load("StarClassLibrary");
+    gSystem->Load("libgsl");
+    gSystem->Load("libgslcblas");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libStDb_Tables.so");
+    gSystem->Load("libgen_Tables.so");
+    gSystem->Load("libgeometry_Tables.so");
+    gSystem->Load("libsim_Tables.so");
+    gSystem->Load("libStarMagField.so");
+    gSystem->Load("libSt_db_Maker.so");
+    gSystem->Load("libSt_g2t.so");
+    gSystem->Load("libSt_geant_Maker.so");
+    gSystem->Load("StarGeneratorUtil.so");
+    gSystem->Load("StarGeneratorEvent.so");
+    gSystem->Load("StarGeneratorBase.so");
+    gSystem->Load("StFemtoDstMakerSph");
+
+    std::cout << "Libraries have been successfully loaded" << std::endl;
+
+    //
+    // Create chain
+    //
+    StChain* mChain = new StChain("StChain");
+    mChain->SetDebug(0);
+    StMuDebug::setLevel(0);
+
+    //
+    // Set StHbtMaker and StHbtManager
+    //
+    StHbtMaker *mHbtMaker = new StHbtMaker();
+    mHbtMaker->SetDebug(0);
+    StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+    //
+    // Set StHbtReader
+    //
+    StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+    Long64_t mNEvents = femtoReader->GetNEvents();
+    mHbtManager->SetEventReader(femtoReader);
+
+    //
+    // Set basic event cut
+    //
+    gregsEventCut *mBasicEventCut = new gregsEventCut();
+    mBasicEventCut->SetCollisionType(0);     //0-pp, 1-AuAu
+    mBasicEventCut->SetEventMult(0, 45);
+    mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+    mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+    mBasicEventCut->SetVertRPos(mPrimVertR);
+
+    //
+    // Set basic kaon single-particle cut
+    //
+    gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+    mBasicTrackCut->SetCharge(mPosCharge);
+    mBasicTrackCut->SetMass(mChargedKaonMass);
+    mBasicTrackCut->SetType(mTrackType);
+    mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+    mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+    mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+    mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+    mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+    mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+    mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+    mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+    mBasicTrackCut->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+    mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+    mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+    mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+    mBasicTrackCut->SetTofMassSqr(mKaonTofMsqr[0],mKaonTofMsqr[1]);
+    mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+    mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+    mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+    mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+    mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+    mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+    mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+    //
+    // Set basic kaon pair cut
+    //
+    gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+    mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+    mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+    mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+    mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+    mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+    Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix, mCurrentPid;
+
+    //
+    // The multiplicity analysis matrix is CHARGE_BINS x MULT_BINS, where
+    // the first element corresponds to the particle charge: (pos, neg),
+    // and the second one corresponds to the RefMult intervals as listed below:
+    // Mult: 0-45, 0-2, 3-4, 5-6, 7-8, 9-10, 11-45
+    //
+    // StHbtAnalysis, cut, CF and monitor declarations
+    //
+    StHbtVertexMultAnalysis *mKaonMultAnalysis[CHARGE_BINS][SPH_BINS];
+
+    //
+    // Cuts
+    //
+    gregsEventCut       *mKaonEventCut[CHARGE_BINS][SPH_BINS];
+    gregsSimpleTrackCut *mKaonTrackCut[CHARGE_BINS][SPH_BINS];
+    gregsTrackPairCut   *mKaonPairCut[CHARGE_BINS][SPH_BINS];
+
+    //
+    // CorrFctns
+    //
+    QinvCorrFctnKt                *mKaonQinvMixKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnOpposHemisphereKt *mKaonQinvOppKtCF[CHARGE_BINS][SPH_BINS];
+    QinvCorrFctnRotationKt        *mKaonQinvRotKtCF[CHARGE_BINS][SPH_BINS];
+
+#ifdef MONITORS
+    //
+    // Cut monitors
+    //
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorPass[CHARGE_BINS][SPH_BINS];
+    fxtEventCutMonitor *mFxtKaonEventCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtTrackCutMonitor *mFxtKaonTrackCutMonitorFail[CHARGE_BINS][SPH_BINS];
+    fxtPairCutMonitor  *mFxtKaonPairCutMonitorFail[CHARGE_BINS][SPH_BINS];
+#endif
+
+    //
+    // Analysis creation
+    //
+    float sphLo, sphHi;
+    mMultLow = 0;
+    mMultHi = 50;
+    mMultBins = 10;
+    mEvents2Mix = 10;
+    for(Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+    {
+        //
+        // sphLo < sphericity <= sphHi
+        //
+        switch (iSph)
+        {
+            case 0:  sphLo = -0.1;     sphHi = 0.1;  break;
+            default: sphLo = iSph*0.1; sphHi = (iSph + 1)*0.1;
+        }
+
+        for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+        {
+            switch(iCharge)
+            {
+                case 0:   mLocalCharge = mPosCharge;    break;
+                default:  mLocalCharge = mNegCharge;    break;
+            };
+
+            //
+            // Create event cut
+            //
+            mKaonEventCut[iCharge][iSph] = new gregsEventCut(*mBasicEventCut);
+            mKaonEventCut[iCharge][iSph]->SetEventMult(mMultLow, mMultHi);
+            mKaonEventCut[iCharge][iSph]->SetSphericityCut(sphLo, sphHi);
+#ifdef MONITORS
+            mFxtKaonEventCutMonitorPass[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonPass", Form("_%d_%d", iCharge, iSph));
+            mFxtKaonEventCutMonitorFail[iCharge][iSph] =
+                new fxtEventCutMonitor("eventKaonFail", Form("_%d_%d", iCharge, iSph));
+            mKaonEventCut[iCharge][iSph]->AddCutMonitor(mFxtKaonEventCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonEventCutMonitorFail[iCharge][iSph]);
+#endif
+
+            //
+            //  Like-sign kaon analysis  
+            //
+            // create analysis
+            mKaonMultAnalysis[iCharge][iSph] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+                                                                           mPrimVertZ[0],
+                                                                           mPrimVertZ[1],
+                                                                           mMultBins,
+                                                                           mMultLow,
+                                                                           mMultHi);
+
+            mKaonMultAnalysis[iCharge][iSph]->SetEventCut(mKaonEventCut[iCharge][iSph]);
+
+            //
+            // Create single-track cuts and add them to the analysis
+            //
+            mKaonTrackCut[iCharge][iSph] = new gregsSimpleTrackCut(*mBasicTrackCut);
+            mKaonTrackCut[iCharge][iSph]->SetCharge(mLocalCharge);
+
+#ifdef MONITORS
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonPass_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph] =
+                new fxtTrackCutMonitor(Form("_trackKaonFail_%d_%d_", iCharge, iSph), mChargedKaonMass);
+            mKaonTrackCut[iCharge][iSph]->AddCutMonitor(mFxtKaonTrackCutMonitorPass[iCharge][iSph],
+                                                        mFxtKaonTrackCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetFirstParticleCut(mKaonTrackCut[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->SetSecondParticleCut(mKaonTrackCut[iCharge][iSph]);
+
+            //
+            // Create pair cut
+            //
+            mKaonPairCut[iCharge][iSph] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+            mFxtKaonPairCutMonitorPass[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonPass_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mFxtKaonPairCutMonitorFail[iCharge][iSph] =
+                new fxtPairCutMonitor(Form("_pairKaonFail_%d_%d_", iCharge, iSph));
+            mFxtKaonPairCutMonitorFail[iCharge][iSph]->SetParticleMass(mChargedKaonMass);
+            mKaonPairCut[iCharge][iSph]->AddCutMonitor(mFxtKaonPairCutMonitorPass[iCharge][iSph],
+                                                       mFxtKaonPairCutMonitorFail[iCharge][iSph]);
+#endif
+            mKaonMultAnalysis[iCharge][iSph]->SetPairCut(mKaonPairCut[iCharge][iSph]);
+
+            //
+            // Correlation function initialization
+            //
+            mKaonQinvMixKtCF[iCharge][iSph] =
+                new QinvCorrFctnKt(Form("hKaonQinvMixKt_%d_%d", iCharge, iSph),
+                                   mQinvBins, mQinvVal[0], mQinvVal[1],
+                                   mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+            mKaonQinvOppKtCF[iCharge][iSph] =
+                new QinvCorrFctnOpposHemisphereKt(Form("hKaonQinvOppKt_%d_%d", iCharge, iSph),
+                                                  mQinvBins, mQinvVal[0], mQinvVal[1],
+                                                  mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                                  mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+            mKaonQinvRotKtCF[iCharge][iSph] =
+                new QinvCorrFctnRotationKt(Form("hKaonQinvRotKt_%d_%d", iCharge, iSph),
+                                           mQinvBins, mQinvVal[0], mQinvVal[1],
+                                           mPairKtBins, mPairKtVal[0], mPairKtVal[1],
+                                           mChargedKaonMass, mKaonPairCut[iCharge][iSph]);
+
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvMixKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvOppKtCF[iCharge][iSph]);
+            mKaonMultAnalysis[iCharge][iSph]->AddCorrFctn(mKaonQinvRotKtCF[iCharge][iSph]);
+
+            mKaonMultAnalysis[iCharge][iSph]->SetNumEventsToMix(mEvents2Mix); // set number of events to mix
+            mHbtManager->AddAnalysis(mKaonMultAnalysis[iCharge][iSph]); // add analysis to the manager
+        }
+    }
+
+    mChain->Init(); // StChain initialization
+
+    //
+    // Start processing
+    //
+    Int_t iRet = 0;
+    Int_t mNEventsProcessed = 0;
+    Float_t mPercentCounter = 0.05;
+    Float_t mProgress = 0.;
+
+    unsigned int mSubCounter = 1;              // time
+    time_t mStartTime, mStopTime, mDiffTime;   // time
+    float mFrac;                               // time
+    mStartTime = time(0);                      // time
+
+    // mNEvents = 15;
+    std::cout << "Total Events: " << mNEvents << std::endl;
+    for(Int_t iEvent=0; iEvent < mNEvents; iEvent++)
+    {
+
+        mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+        mNEventsProcessed++;
+
+        if(mProgress >= mPercentCounter)
+        {
+            mPercentCounter += 0.05;
+            mStopTime = time(0);
+            mDiffTime = difftime(mStopTime, mStartTime);
+            mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+            mSubCounter++;
+            std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+                << std::endl;
+        }
+
+        mChain->Clear();
+        iRet = mChain->Make(iEvent);
+        if(iRet!=0)
+        {
+            std::cout << "[ERROR] - iRet: " << iRet << std::endl;
+            break;
+        }
+    }
+
+    mChain->Finish();
+
+    TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+
+    //
+    // Write hbt histograms to the file
+    //
+    std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+    for(Int_t iCharge = 0; iCharge < CHARGE_BINS; iCharge++)
+    {
+        for (Int_t iSph = 0; iSph < SPH_BINS; iSph++)
+        {
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexZ()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->RefMult()->Write();
+            mFxtKaonEventCutMonitorPass[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexZ()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VertexYvsVertexX()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->RefMult()->Write();
+            mFxtKaonEventCutMonitorFail[iCharge][iSph]->VpdVzDiff()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Nhits()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->P()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Pt()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsDedx()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->Rapidity()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtKaonTrackCutMonitorPass[iCharge][iSph]->PvsInvBeta()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->DCAGlobal()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Nhits()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->P()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Pt()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaPion()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaKaon()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PtVsNsigmaProton()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsDedx()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->Rapidity()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PseudoRapidity()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsMassSqr()->Write();
+            mFxtKaonTrackCutMonitorFail[iCharge][iSph]->PvsInvBeta()->Write();
+
+            mFxtKaonPairCutMonitorPass[iCharge][iSph]->Kt()->Write();