Browse Source

First commit

Parfenov_Petr_Evgen'evich 6 months ago
commit
b57e5eefb6
12 changed files with 1369 additions and 0 deletions
  1. 83 0
      README.md
  2. 62 0
      macro/CMakeLists.txt
  3. 271 0
      macro/CalcPID.C
  4. 64 0
      macro/Constants.h
  5. 287 0
      macro/FemtoDstAnalyzer_PID.C
  6. 398 0
      macro/FemtoDstAnalyzer_PID_QA.C
  7. 60 0
      macro/mainPID.cpp
  8. 70 0
      macro/mainPID_QA.cpp
  9. 36 0
      scripts/GenerateLists.sh
  10. 12 0
      scripts/run.sh
  11. 22 0
      scripts/start.sh
  12. 4 0
      set_env.sh

+ 83 - 0
README.md

@@ -0,0 +1,83 @@
+# Quick information
+***
+
+It is a simple example of scripts/macros sets suitable for work on the MEPhI (basov) cluster.
+
+## Setting environment
+
+Copy from a git repo:
+
+        git clone https://devel.mephi.ru/PEParfenov/hpc_scripts.git
+
+Change paths accordingly:
+
+* 1 In `set_env.sh`: Change `ST_FEMTO_DST_INC_DIR` to standard one. It's the directory where `libStFemtoDst.so` is stored.
+* 2 In `scripts/run.sh`: change path to a last line to one where your build will be.
+* 3 In `scripts/start.sh`: change path to output directory (`OUTPUT_DIR`).
+
+## Installation
+
+Source required environment variables:
+
+        cd pid_scripts/
+        . set_env.sh
+
+Make new build directory. For example:
+
+        mkdir build/
+        cd build/
+
+Generate makefile & install:
+
+        cmake ../macro/
+        make
+
+Keep in mind that one has to do `make` command after changing `FemtoDstAnalyzer.C` in order to compile new code.
+
+## Generate filelists
+
+Use `GenerateLists.sh` to make filelists:
+
+        . GenerateLists.sh FEMTODST_DIR N_FILES_IN_LIST
+
+where `FEMTODST_DIR` - path to the directory with `femtoDst.root` files.
+And `N_FILES_IN_LIST` denotes the maximum number of `femtoDst.root` files in each filelist.
+Basic example:
+
+          . GenerateLists.sh /mnt/pool/rhic/2/nigmatkulov/femtoDst/auau/200gev/12135/ 100
+
+Resulting filelists will be in the `hpc_scripts/lists/` directory.
+
+## Usage
+
+### Interactive mode
+
+To use `FemtoDstAnalyzer_PID.C` in interactive mode:
+
+        cd build/
+        ./pid -i INPUTFILE -o OUTPUTFILE
+
+To use `FemtoDstAnalyzer_PID.C` in interactive mode:
+
+        cd build/
+        ./pid -i INPUTFILE -o OUTPUTFILE -pid PID_FILE
+
+where `INPUTFILE` - is input file or filelist with `femtoDst.root`.
+`OUTPUTFILE` - is resulting root file and `PID_FILE` - is fitted sigma values in m^{2} distribution.
+Basic workflow example:
+
+        ./pid -i ../lists/StRuns1.list -o ./test_pid.root
+	root.exe '../macro/CalcPID.C("test_pid.root","test_pid_fit.root",1)'
+	./pidQA -i ../lists/StRuns1.list -o ./test_pid_QA.root -pid ./test_pid_fit.root
+
+### Batch mode
+
+To send jobs to basov cluster, use `scripts/start.sh`:
+
+        . start.sh INPUT_FILELIST_DIR MODE_NUM
+
+where `INPUT_FILELIST_DIR` - is the directory where you store filelists.
+`MOD_NUM` - is mode number: 0 - run `pid`, 1 - run `pidQA`.
+Basic example:
+
+        . start.sh /mnt/pool/rhic/4/parfenovpeter/STAR/Analysis/pid_scripts/lists 0

+ 62 - 0
macro/CMakeLists.txt

@@ -0,0 +1,62 @@
+# CMakeLists.txt for FAS package. It creates a library with dictionary and a main program
+cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
+project(FemtoDstAnalyzer)
+
+# You need to tell CMake where to find the ROOT installation. This can be done
+# in a number of ways:
+#   - ROOT built with classic configure/make use the provided
+#   $ROOTSYS/etc/cmake/FindROOT.cmake
+#   - ROOT built with CMake. Add in CMAKE_PREFIX_PATH the installation prefix
+#   for ROOT
+if (ROOT_CMAKE)
+  list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS})
+else (ROOT_CMAKE)
+  set(ROOT_PREINSTALLED "/usr/lib64/Geant4-10.3.0/Modules")
+  list(APPEND CMAKE_MODULE_PATH ${ROOT_PREINSTALLED})
+endif (ROOT_CMAKE)
+
+#---Locate the ROOT package and defines a number of variables (e.g. ROOT_INCLUDE_DIRS)
+if (ROOT_CMAKE)
+  find_package(ROOT REQUIRED COMPONENTS MathCore RIO Hist Tree Net
+	  HINTS ENV{PATH})
+else (ROOT_CMAKE)
+  find_package(ROOT REQUIRED COMPONENTS MathCore RIO Hist Tree Net)
+endif (ROOT_CMAKE)
+
+#---Define useful ROOT functions and macros (e.g. ROOT_GENERATE_DICTIONARY)
+if (ROOT_CMAKE)
+  include(${ROOT_USE_FILE})
+endif (ROOT_CMAKE)
+
+
+#---Define _VANILLA_ROOT_ variable needed for the project
+#add_definitions(-D_VANILLA_ROOT_)
+
+set(ST_FEMTO_DST_LIB libStFemtoDst.so)
+set(ST_FEMTO_DST_CUSTOM_DIR $ENV{ST_FEMTO_DST_BUILD_DIR})
+set(ST_FEMTO_DST_CUSTOM_INC_DIR $ENV{ST_FEMTO_DST_INC_DIR})
+
+find_library(ST_FEMTO_DST REQUIRED
+	NAMES ${ST_FEMTO_DST_LIB}
+	HINTS ${ST_FEMTO_DST_CUSTOM_DIR})
+
+set(INCLUDE_DIRECTORIES
+  ${CMAKE_SOURCE_DIR}
+  ${ST_FEMTO_DST_CUSTOM_INC_DIR}
+  ${ROOT_INCLUDE_DIRS}
+)
+
+include_directories(${INCLUDE_DIRECTORIES})
+set (CMAKE_CXX_STANDARD 11)
+add_definitions(${ROOT_CXX_FLAGS})
+
+set(CMAKE_BUILD_TYPE Debug)
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O2")
+
+#---Make executable
+
+add_executable(pid mainPID.cpp)
+target_link_libraries(pid ${ST_FEMTO_DST} ${ROOT_LIBRARIES})
+
+add_executable(pidQA mainPID_QA.cpp)
+target_link_libraries(pidQA ${ST_FEMTO_DST} ${ROOT_LIBRARIES})

+ 271 - 0
macro/CalcPID.C

@@ -0,0 +1,271 @@
+#include <iostream>
+#include <TFile.h>
+#include <TF2.h>
+#include <TH2F.h>
+#include <TMath.h>
+
+const int NptBins = 24;
+
+const double pionMassSqr = TMath::Power(0.13957061, 2);
+const double kaonMassSqr = TMath::Power(0.493677, 2);
+const double protMassSqr = TMath::Power(0.938272081, 2);
+
+void Fit2Dgaus(const Char_t *inFileName, const Char_t *outFileName)
+{
+
+  const double ptBin[NptBins + 1] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8, 4., 4.5, 5., 5.5, 6.};
+
+  const double pionNsigMean[NptBins] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
+  const double kaonNsigMean[NptBins] = {12., 6., 3., 1., -0.3, -1., -1.2, -1.5, -1.5, -1.5, -1.7, -2., -2., -2., -2., -2., -2., -2., -2., -2., -2., -2., -2., -2.};
+  const double protNsigMean[NptBins] = {22., 15., 9., 5., 3., 1.5, 0.2, -0.5, -1., -1.5, -1.8, -2., -2.2, -2.2, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5};
+
+  const double pionNsigWidth[NptBins] = {4., 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
+  const double kaonNsigWidth[NptBins] = {8., 5., 4., 4., 3.7, 3., 3., 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
+  const double protNsigWidth[NptBins] = {8., 10., 7., 5., 4., 3.5, 3., 2.6, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
+
+  const std::pair<double,double> pionMsqrRange [NptBins] = {{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},
+                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.1,0.14},
+                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.036,0.019},
+                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},
+                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},
+                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.039,0.078}};
+
+  const std::pair<double,double> kaonMsqrRange [NptBins] = {{0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
+                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.13,0.4},
+                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
+                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
+                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
+                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38}};
+
+  const std::pair<double,double> protMsqrRange [NptBins] = {{0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
+                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
+                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
+                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
+                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
+                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4}};
+
+  const std::pair<double,double> tripleNsigRange [NptBins] = {{-6., 28.},{-6., 28.},{-6., 28.},{-6., 28.},
+                                                              {-6., 28.},{-6., 28.},{-6., 28.},{-5., 4.},
+                                                              {-5., 4.},{-6., 28.},{-6., 28.},{-5., 3.},
+                                                              {-5., 3.},{-5., 3.},{-5., 3.},{-5., 3.},
+                                                              {-5., 3.},{-5., 3.},{-5., 3.},{-5., 3.},
+                                                              {-5., 2.},{-4.5, 2.},{-4., 2.},{-4., 2.}};
+
+  const std::pair<double,double> tripleMsqrRange [NptBins] = {{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
+                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.4,1.1},
+                                                              {-0.3,1.3},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
+                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
+                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
+                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5}};
+
+  TFile *fi = new TFile(inFileName, "read");
+
+  TH2F *hNsigMsqr[NptBins];
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigMsqr[i] = (TH2F *)fi->Get(Form("hNsigPionMSqr%i", i));
+  }
+
+  TF2 *gausPion[NptBins];
+  TF2 *gausKaon[NptBins];
+  TF2 *gausProton[NptBins];
+  TF2 *gausTriple[NptBins];
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    gausPion[i] = new TF2(Form("gausPion%i", i), "xygaus", pionNsigMean[i] - pionNsigWidth[i], pionNsigMean[i] + pionNsigWidth[i], pionMsqrRange[i].first, pionMsqrRange[i].second);
+    gausKaon[i] = new TF2(Form("gausKaon%i", i), "xygaus", kaonNsigMean[i] - kaonNsigWidth[i], kaonNsigMean[i] + kaonNsigWidth[i], kaonMsqrRange[i].first, kaonMsqrRange[i].second);
+    gausProton[i] = new TF2(Form("gausProton%i", i), "xygaus", protNsigMean[i] - protNsigWidth[i], protNsigMean[i] + protNsigWidth[i], protMsqrRange[i].first, protMsqrRange[i].second);
+    gausTriple[i] = new TF2(Form("gausTriple%i", i), "xygaus(0)+xygaus(5)+xygaus(10)", tripleNsigRange[i].first, tripleNsigRange[i].second, tripleMsqrRange[i].first, tripleMsqrRange[i].second);
+  }
+
+  // gausPion[0] = new TF2(Form("gausPion%i",0),"xygaus",-3.,3.,0.,0.4);
+  // gausKaon[0] = new TF2(Form("gausKaon%i",0),"xygaus",4.,20.,0.2,0.28);
+  // gausProton[0] = new TF2(Form("gausProton%i",0),"xygaus",6.8,32.85,0.54,1.25);
+  // gausProton[0]->SetParameters(7.97,20.85,0.888,0.04366);
+  // gausTriple[0] = new TF2(Form("gausTriple%i",0),"xygaus(0)+xygaus(5)+xygaus(10)",-5.,32.85,-0.5,2.);
+
+  int TestStop = 4;
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigMsqr[i]->Fit(gausPion[i], "R0");
+    hNsigMsqr[i]->Fit(gausKaon[i], "R0");
+    hNsigMsqr[i]->Fit(gausProton[i], "R0");
+
+    Double_t par[15];
+    gausPion[i]->GetParameters(&par[0]);
+    gausKaon[i]->GetParameters(&par[5]);
+    gausProton[i]->GetParameters(&par[10]);
+    gausTriple[i]->SetParameters(par);
+
+    std::cout << "Fit plot: " << i << std::endl;
+    hNsigMsqr[i]->Fit(gausTriple[i], "R");
+  }
+
+  TFile *fo = new TFile(outFileName, "recreate");
+  fo->cd();
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigMsqr[i]->Write();
+  }
+  for (int i = 0; i < NptBins; i++)
+  {
+    gausPion[i]->Write();
+    gausKaon[i]->Write();
+    gausProton[i]->Write();
+  }
+  for (int i = 0; i < NptBins; i++)
+  {
+    gausTriple[i]->Write();
+  }
+}
+
+void Fit1DgausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
+{
+  TFile *fi = new TFile(inFileName, "read");
+
+  TH2F *hNsigMsqrPion[NptBins];
+  TH2F *hNsigMsqrKaon[NptBins];
+  TH2F *hNsigMsqrProton[NptBins];
+
+  TH1F *hMsqrPion[NptBins];
+  TH1F *hMsqrKaon[NptBins];
+  TH1F *hMsqrProton[NptBins];
+
+  TF1 *gausPion[NptBins];
+  TF1 *gausKaon[NptBins];
+  TF1 *gausProton[NptBins];
+  
+  TF1 *gausTriplePion[NptBins];
+  TF1 *gausTripleKaon[NptBins];
+  TF1 *gausTripleProton[NptBins];
+
+  TH1F *hMsqrSigmaPion = new TH1F("hMsqrSigmaPion","hMsqrSigmaPion",NptBins,0.15,6.15);
+  TH1F *hMsqrSigmaKaon = new TH1F("hMsqrSigmaKaon","hMsqrSigmaKaon",NptBins,0.15,6.15);
+  TH1F *hMsqrSigmaProton = new TH1F("hMsqrSigmaProton","hMsqrSigmaProton",NptBins,0.15,6.15);
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigMsqrPion[i] = (TH2F *)fi->Get(Form("hNsigPionMSqr%i", i));
+    hNsigMsqrKaon[i] = (TH2F *)fi->Get(Form("hNsigKaonMSqr%i", i));
+    hNsigMsqrProton[i] = (TH2F *)fi->Get(Form("hNsigProtonMSqr%i", i));
+
+    gausPion[i] = new TF1(Form("gausPion%i",i),"gaus",pionMassSqr*0.,pionMassSqr*2.);
+    gausKaon[i] = new TF1(Form("gausKaon%i",i),"gaus",kaonMassSqr*0.5,kaonMassSqr*1.5);
+    gausProton[i] = new TF1(Form("gausProton%i",i),"gaus",protMassSqr*0.5,protMassSqr*1.5);
+    
+    gausTriplePion[i] = new TF1(Form("gausTriplePion%i",i),"gaus(0)+gaus(3)+gaus(6)",-0.5,1.5);
+    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i",i),"gaus(0)+gaus(3)+gaus(6)",-0.5,1.5);
+    gausTripleProton[i] = new TF1(Form("gausTripleProton%i",i),"gaus(0)+gaus(3)+gaus(6)",-0.5,1.5);
+  }
+
+  int bin1 = 0, bin2 = 0;
+  Double_t parPion[9], parKaon[9], parProton[9];
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    bin1 = 651; //hNsigMsqrPion[i]->FindBin(-2.5);
+    bin2 = 750; //hNsigMsqrPion[i]->FindBin(2.5);
+    hMsqrPion[i] = (TH1F*) hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i",i),bin1,bin2);
+    hMsqrKaon[i] = (TH1F*) hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i",i),bin1,bin2);
+    hMsqrProton[i] = (TH1F*) hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i",i),bin1,bin2);
+  }
+
+  //Pion
+  for (int i = 0; i < NptBins; i++)
+  {
+    hMsqrPion[i]->Fit(gausPion[i],"RQ0");
+    hMsqrPion[i]->Fit(gausKaon[i],"RQ0");
+    hMsqrPion[i]->Fit(gausProton[i],"RQ0");
+
+    gausPion[i]->GetParameters(&parPion[0]);
+    gausKaon[i]->GetParameters(&parPion[3]);
+    gausProton[i]->GetParameters(&parPion[6]);
+
+    gausTriplePion[i]->SetParameters(parPion);
+
+    hMsqrPion[i]->Fit(gausTriplePion[i],"R");
+    hMsqrSigmaPion->SetBinContent(i,gausTriplePion[i]->GetParameter(2));
+    hMsqrSigmaPion->SetBinError(i,gausTriplePion[i]->GetParError(2));
+  }
+
+  //Kaon
+  for (int i = 0; i < NptBins; i++)
+  {
+    hMsqrKaon[i]->Fit(gausPion[i],"RQ0");
+    hMsqrKaon[i]->Fit(gausKaon[i],"RQ0");
+    hMsqrKaon[i]->Fit(gausProton[i],"RQ0");
+
+    gausPion[i]->GetParameters(&parKaon[0]);
+    gausKaon[i]->GetParameters(&parKaon[3]);
+    gausProton[i]->GetParameters(&parKaon[6]);
+
+    gausTripleKaon[i]->SetParameters(parKaon);
+
+    hMsqrKaon[i]->Fit(gausTripleKaon[i],"R");
+    hMsqrSigmaKaon->SetBinContent(i,gausTripleKaon[i]->GetParameter(5));
+    hMsqrSigmaKaon->SetBinError(i,gausTripleKaon[i]->GetParError(5));
+  }
+
+  //Proton
+  for (int i = 0; i < NptBins; i++)
+  {
+    hMsqrProton[i]->Fit(gausPion[i],"RQ0");
+    hMsqrProton[i]->Fit(gausKaon[i],"RQ0");
+    hMsqrProton[i]->Fit(gausProton[i],"RQ0");
+
+    gausPion[i]->GetParameters(&parProton[0]);
+    gausKaon[i]->GetParameters(&parProton[3]);
+    gausProton[i]->GetParameters(&parProton[6]);
+
+    gausTripleProton[i]->SetParameters(parProton);
+
+    hMsqrProton[i]->Fit(gausTripleProton[i],"R");
+    hMsqrSigmaProton->SetBinContent(i,gausTripleProton[i]->GetParameter(8));
+    hMsqrSigmaProton->SetBinError(i,gausTripleProton[i]->GetParError(8));
+  }
+
+  TF1 *polPion = new TF1("polPion","pol8",0.15,4.9);
+  TF1 *polKaon = new TF1("polKaon","pol8",0.15,4.9);
+  TF1 *polProton = new TF1("polProton","pol8",0.15,4.9);
+
+  hMsqrSigmaPion->Fit(polPion,"R");
+  hMsqrSigmaKaon->Fit(polKaon,"R");
+  hMsqrSigmaProton->Fit(polProton,"R");
+
+  TFile *fo = new TFile(outFileName, "recreate");
+  fo->cd();
+  // for (int i = 0; i < NptBins; i++)
+  // {
+  //   hMsqrPion[i]->Write();
+  //   hMsqrKaon[i]->Write();
+  //   hMsqrProton[i]->Write();
+  // }
+  // for (int i = 0; i < NptBins; i++)
+  // {
+  //   gausTriplePion[i]->Write();
+  //   gausTripleKaon[i]->Write();
+  //   gausTripleProton[i]->Write();
+  // }
+  hMsqrSigmaPion->Write();
+  hMsqrSigmaKaon->Write();
+  hMsqrSigmaProton->Write();
+
+  polPion->Write();
+  polKaon->Write();
+  polProton->Write();
+}
+
+void CalcPID(const Char_t *inFileName, const Char_t *outFileName, const Int_t mode)
+{
+  if (mode==0)
+  {
+    Fit2Dgaus(inFileName, outFileName);
+  }
+  if (mode==1)
+  {
+    Fit1DgausMassSqr(inFileName, outFileName);
+  }
+}

+ 64 - 0
macro/Constants.h

@@ -0,0 +1,64 @@
+#include <TMath.h>
+/////////////////////////////////////CONSTANTS/////////////////////////////////////////////////
+
+//Used runIds
+const int runId_min = 12126101;
+const int runId_max = 12134068;
+
+// Constants
+const Double_t cutVtxZ = 40.;
+const std::map<Float_t, Double_t> cutVtxZEnergy = {{7.7, 70.}, {11.5, 50.}, {19.6, 70}, {27., 70.}, {39., 40.}, {62., 40.}, {200., 30.}};
+const int fArraySize = 2;  // vertex pos/neg
+const int fNharmonics = 2; // v2 and v3
+const int fNBins = runId_max - runId_min + 1;
+
+//Event cuts
+const Double_t cutVtxR = 2.;
+const Double_t cutVtxZvpd = 3.;
+const Int_t cutNTofPoints = 2.;
+const Double_t cutVpdVz = 3.;
+
+//Track cuts
+const std::map<Float_t, Double_t> cutDCA = {{7.7, 1.}, {11.5, 1.}, {19.6, 1.}, {27., 1.}, {39., 1.}, {62., 1.}, {200., 3.}};
+const Double_t cutDCA_PID = 1.;
+const Double_t cutEta = 1.;
+const Int_t    cutNhits = 15;
+const Double_t cutPtotMin = 0.1;
+const Double_t cutNhitsPoss = 0.;
+const Double_t cutNhitsRatio = 0.51;
+const std::map<Float_t, Double_t> cutPtMin = {{7.7, 0.2}, {11.5, 0.2}, {19.6, 0.2}, {27., 0.2}, {39., 0.2}, {62., 0.2}, {200., 0.15}};
+const Double_t cutPtMax = 2.;
+const Double_t cutPMax = 10.;
+
+//PID cuts
+const Double_t cutMass2Min = -10.;
+const Double_t cutNsigPID = 2.5;
+
+//ReCentering
+const Double_t cutPtWeightEP = 2.;
+const std::vector<Double_t> cutEtaGap = {0.05,0.1,0.2,0.5};
+const Int_t    NEtaGaps = cutEtaGap.size();
+
+//Shift Corrections
+const Int_t NShiftOrderMax = 4;
+const Int_t cutNcountQvFull = 4;
+const Int_t cutNcountQvGap = 2;
+const Int_t cutNcountQvFullEast = 0;
+const Int_t cutNcountQvFullWest = 0;
+const std::vector<Double_t> shiftOrderCoeff2 = {2.0, 4.0, 6.0, 8.0, 10.0};
+const std::vector<Double_t> shiftOrderCoeff3 = {3.0, 6.0, 9.0, 12.0, 15.0};
+
+//Energy of the collision
+const Double_t energy = 200.;
+
+//PID
+const double pionMassSqr = TMath::Power(0.13957061, 2);
+const double kaonMassSqr = TMath::Power(0.493677, 2);
+const double protMassSqr = TMath::Power(0.938272081, 2);
+
+const double NpidFactor = 2.5;
+
+const double Nquarks[3] = {2.,2.,3.}; // pion, kaon, proton
+const double MassSqr[3] = {pionMassSqr, kaonMassSqr, protMassSqr};
+
+/////////////////////////////////////CONSTANTS/////////////////////////////////////////////////

+ 287 - 0
macro/FemtoDstAnalyzer_PID.C

@@ -0,0 +1,287 @@
+/**
+ * \brief Example of how to read a file (list of files) using StFemtoEvent classes
+ *
+ * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
+ * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
+ * inFile.list) as an input, and preform physics analysis
+ *
+ * \author Grigory Nigmatkulov
+ * \date May 29, 2018
+ */
+
+// This is needed for calling standalone classes
+#define _VANILLA_ROOT_
+
+// C++ headers
+#include <iostream>
+#include <vector>
+#include <map>
+
+// ROOT headers
+#include "TROOT.h"
+#include "TFile.h"
+#include "TChain.h"
+#include "TVector2.h"
+#include "TVector3.h"
+#include "TTree.h"
+#include "TSystem.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TMath.h"
+#include "TProfile2D.h"
+#include "TStopwatch.h"
+
+// FemtoDst headers
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDstReader.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDst.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoEvent.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoTrack.h"
+
+// Load libraries (for ROOT_VERSTION_CODE >= 393215)
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
+R__LOAD_LIBRARY(/mnt/pool/rhic/4/parfenovpeter/STAR/build/libStFemtoDst.so)
+#endif
+
+#include "Constants.h"
+
+// inFile - is a name of name.FemtoDst.root file or a name
+//          of a name.lis(t) files, that contains a list of
+//          name1.FemtoDst.root, name2.FemtoDst.root, ... files
+
+//Used functions (see them below)
+Bool_t isGoodEvent(StFemtoEvent *const &event);
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Double_t GetWeight(StFemtoTrack *const &track);
+TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
+Double_t AngleShift(Double_t Psi, Double_t order);
+
+//_________________
+void FemtoDstAnalyzer_PID(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
+                                         const Char_t *outFile = "./oPIDTest.root")
+{
+
+  std::cout << "Hi! Lets do some physics, Master!" << std::endl;
+  TStopwatch timer;
+  timer.Start();
+
+#if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
+  gSystem->Load("../libStFemtoDst.so");
+#endif
+
+  StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
+  femtoReader->Init();
+
+  // This is a way if you want to spead up IO
+  std::cout << "Explicit read status for some branches" << std::endl;
+  femtoReader->SetStatus("*", 0);
+  femtoReader->SetStatus("Event", 1);
+  femtoReader->SetStatus("Track", 1);
+  std::cout << "Status has been set" << std::endl;
+
+  std::cout << "Now I know what to read, Master!" << std::endl;
+
+  if (!femtoReader->chain())
+  {
+    std::cout << "No chain has been found." << std::endl;
+  }
+  Long64_t eventsInTree = femtoReader->tree()->GetEntries();
+  std::cout << "eventsInTree: " << eventsInTree << std::endl;
+  Long64_t events2read = femtoReader->chain()->GetEntries();
+
+  std::cout << "Number of events to read: " << events2read << std::endl;
+
+
+  const int NptBins = 24;
+  const double ptBin [NptBins+1] = {0.,0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.2,3.4,3.6,3.8,4.,4.5,5.,5.5,6.};
+  TH2F *hNsigPionMSqr[NptBins];
+  TH2F *hNsigKaonMSqr[NptBins];
+  TH2F *hNsigProtonMSqr[NptBins];
+
+  //Initialization
+  for (int i=0; i<NptBins; i++)
+  {
+    hNsigPionMSqr[i] = new TH2F(Form("hNsigPionMSqr%i",i) ,Form("n#sigma_{#pi} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{#pi};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
+    hNsigKaonMSqr[i] = new TH2F(Form("hNsigKaonMSqr%i",i) ,Form("n#sigma_{K} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{K};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
+    hNsigProtonMSqr[i] = new TH2F(Form("hNsigProtonMSqr%i",i) ,Form("n#sigma_{p} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{p};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
+  }
+
+  // Loop over events
+  for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
+  {
+
+    if ((iEvent+1) % 1000 == 0){
+      std::cout << "Working on event #[" << (iEvent + 1)
+                << "/" << events2read << "]" << std::endl;
+    }
+
+    Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
+    if (!readEvent)
+    {
+      std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
+      break;
+    }
+
+    // Retrieve femtoDst
+    StFemtoDst *dst = femtoReader->femtoDst();
+
+    // Retrieve event information
+    StFemtoEvent *event = dst->event();
+    if (!event)
+    {
+      std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
+      break;
+    }
+
+    TVector3 pVtx = event->primaryVertex();
+
+    // Event selection
+    if ( !isGoodEvent(event) ) continue;
+
+    // Track analysis
+    Int_t nTracks = dst->numberOfTracks();
+    
+
+    // Track loop
+    for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
+    {
+
+      // Retrieve i-th femto track
+      StFemtoTrack *femtoTrack = dst->track(iTrk);
+
+      if (!femtoTrack)
+        continue;
+
+      // Must be a primary track
+      if (!femtoTrack->isPrimary())
+        continue;
+
+      //Track selection
+      if (!isGoodTrackFlow(femtoTrack, energy, pVtx)) continue;
+
+      // Determine pt bin for PID plots
+      int i_pt = -1;
+      for (int i=0;i<NptBins;i++)
+      {
+        if (femtoTrack->pt() >= (ptBin[i] + cutPtMin.at(energy)) && femtoTrack->pt() < (ptBin[i+1] + cutPtMin.at(energy)))
+	{
+	  i_pt = i;
+	}
+      }
+      if (i_pt == -1) continue;
+
+      // Check if track has TOF signal
+      if (femtoTrack->isTofTrack())
+      {
+        if (femtoTrack->gDCA(pVtx).Mag() >= 1.) continue;
+
+        hNsigPionMSqr[i_pt]->Fill(femtoTrack->nSigmaPion(),femtoTrack->massSqr());
+        hNsigKaonMSqr[i_pt]->Fill(femtoTrack->nSigmaKaon(),femtoTrack->massSqr());
+        hNsigProtonMSqr[i_pt]->Fill(femtoTrack->nSigmaProton(),femtoTrack->massSqr());
+
+      } //if( isTofTrack() )
+
+    } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
+
+  } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
+
+  femtoReader->Finish();
+
+  TFile *output = new TFile(outFile,"recreate");
+
+  for (int i=0; i<NptBins; i++)
+  {
+    hNsigPionMSqr[i] -> Write();
+    hNsigKaonMSqr[i] -> Write();
+    hNsigProtonMSqr[i] -> Write();
+  }
+
+  output->Close();
+
+  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
+            << std::endl;
+  timer.Stop();
+  timer.Print();
+}
+
+Bool_t isGoodEvent(StFemtoEvent *const &event)
+{
+  if (!event) return false;
+  if (event == nullptr) return false;
+
+  if (event->primaryVertex().Perp() > cutVtxR) return false;
+  if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy)) return false;
+
+  if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz) return false;
+
+  return true;
+}
+
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track) return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits) return false;
+  if (track->nHitsPoss() <= cutNhitsPoss) return false;
+  if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
+  if (TMath::Abs(track->eta()) >= cutEta) return false;
+  
+  if (track->pt() <= cutPtMin.at(_energy)) return false;
+  if (track->pt() > cutPtMax) return false;
+  if (track->ptot() > cutPMax) return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
+  return true;
+}
+
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track) return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits) return false;
+  if (track->nHitsPoss() <= cutNhitsPoss) return false;
+  if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
+  if (TMath::Abs(track->eta()) >= cutEta) return false;
+  
+  if (track->pt() <= cutPtMin.at(_energy)) return false;
+  //if (track->pt() > cutPtMax) return false;
+  if (track->ptot() > cutPMax) return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
+  return true;
+}
+
+Double_t GetWeight(StFemtoTrack *const &track)
+{
+  Double_t weight;
+  if (track->pt() <= cutPtWeightEP)
+  {
+    weight = track->pt();
+  }
+  else
+  {
+    weight = cutPtWeightEP;
+  }
+  return weight;
+}
+
+TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
+{
+  TVector2 qv(0.,0.);
+
+  qv.Set(TMath::Cos(_harm*track->phi()),TMath::Sin(_harm*track->phi()));
+
+  return qv;
+}
+
+Double_t AngleShift(Double_t Psi, Double_t order)
+{
+  Double_t PsiCorr = Psi;
+  if (PsiCorr > TMath::Pi()/order)
+  {
+    PsiCorr = Psi - 2.*TMath::Pi()/order;
+  }
+  if (PsiCorr < -1.*TMath::Pi()/order)
+  {
+    PsiCorr = Psi + 2.*TMath::Pi()/order;
+  }
+  return PsiCorr;
+}

+ 398 - 0
macro/FemtoDstAnalyzer_PID_QA.C

@@ -0,0 +1,398 @@
+/**
+ * \brief Example of how to read a file (list of files) using StFemtoEvent classes
+ *
+ * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
+ * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
+ * inFile.list) as an input, and preform physics analysis
+ *
+ * \author Grigory Nigmatkulov
+ * \date May 29, 2018
+ */
+
+// This is needed for calling standalone classes
+#define _VANILLA_ROOT_
+
+// C++ headers
+#include <iostream>
+#include <vector>
+#include <map>
+
+// ROOT headers
+#include "TROOT.h"
+#include "TFile.h"
+#include "TChain.h"
+#include "TVector2.h"
+#include "TVector3.h"
+#include "TTree.h"
+#include "TSystem.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH2.h"
+#include "TMath.h"
+#include "TProfile2D.h"
+#include "TStopwatch.h"
+
+// FemtoDst headers
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDstReader.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDst.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoEvent.h"
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoTrack.h"
+
+// Load libraries (for ROOT_VERSTION_CODE >= 393215)
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
+R__LOAD_LIBRARY(/mnt/pool/rhic/4/parfenovpeter/STAR/build/libStFemtoDst.so)
+#endif
+
+#include "Constants.h"
+
+// inFile - is a name of name.FemtoDst.root file or a name
+//          of a name.lis(t) files, that contains a list of
+//          name1.FemtoDst.root, name2.FemtoDst.root, ... files
+
+//Used functions (see them below)
+Bool_t isGoodEvent(StFemtoEvent *const &event);
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Double_t GetWeight(StFemtoTrack *const &track);
+TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
+Double_t AngleShift(Double_t Psi, Double_t order);
+
+//_________________
+void FemtoDstAnalyzer_PID_QA(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
+                            const Char_t *outFile = "./oPIDTest.root",
+                            const Char_t *pidFile = "")
+{
+
+  std::cout << "Hi! Lets do some physics, Master!" << std::endl;
+  TStopwatch timer;
+  timer.Start();
+
+#if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
+  gSystem->Load("../libStFemtoDst.so");
+#endif
+
+  StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
+  femtoReader->Init();
+
+  // This is a way if you want to spead up IO
+  std::cout << "Explicit read status for some branches" << std::endl;
+  femtoReader->SetStatus("*", 0);
+  femtoReader->SetStatus("Event", 1);
+  femtoReader->SetStatus("Track", 1);
+  std::cout << "Status has been set" << std::endl;
+
+  std::cout << "Now I know what to read, Master!" << std::endl;
+
+  if (!femtoReader->chain())
+  {
+    std::cout << "No chain has been found." << std::endl;
+  }
+  Long64_t eventsInTree = femtoReader->tree()->GetEntries();
+  std::cout << "eventsInTree: " << eventsInTree << std::endl;
+  Long64_t events2read = femtoReader->chain()->GetEntries();
+
+  std::cout << "Number of events to read: " << events2read << std::endl;
+
+  //Read PID fit parameters
+  TFile *fiPID = new TFile(pidFile, "read");
+  fiPID->cd();
+  TF1 *sigMsqrPion = (TF1 *)fiPID->Get("polPion");
+  TF1 *sigMsqrKaon = (TF1 *)fiPID->Get("polKaon");
+  TF1 *sigMsqrProton = (TF1 *)fiPID->Get("polProton");
+
+  const int NptBins = 24;
+  const double ptBin[NptBins + 1] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8, 4., 4.5, 5., 5.5, 6.};
+  TH2F *hNsigPionMSqr[NptBins];
+  TH2F *hNsigKaonMSqr[NptBins];
+  TH2F *hNsigProtonMSqr[NptBins];
+
+  TH1F *hNsigPionMSqrAfter[NptBins];
+  TH1F *hNsigKaonMSqrAfter[NptBins];
+  TH1F *hNsigProtonMSqrAfter[NptBins];
+
+  //Initialization
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigPionMSqr[i] = new TH2F(Form("hNsigPionMSqr%i", i), Form("n#sigma_{#pi} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{#pi};m^{2} [GeV/c^{2}]^{2}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -35., 35., 1400, -2., 5.);
+    hNsigKaonMSqr[i] = new TH2F(Form("hNsigKaonMSqr%i", i), Form("n#sigma_{K} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{K};m^{2} [GeV/c^{2}]^{2}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -35., 35., 1400, -2., 5.);
+    hNsigProtonMSqr[i] = new TH2F(Form("hNsigProtonMSqr%i", i), Form("n#sigma_{p} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{p};m^{2} [GeV/c^{2}]^{2}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -35., 35., 1400, -2., 5.);
+
+    hNsigPionMSqrAfter[i] = new TH1F(Form("hNsigPionMSqrAfter%i", i), Form("M^{2} for %.2f < p_{T} < %.2f GeV/c;m^{2} [GeV/c^{2}]^{2};N_{counts}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -2., 5.);
+    hNsigKaonMSqrAfter[i] = new TH1F(Form("hNsigKaonMSqrAfter%i", i), Form("M^{2} for %.2f < p_{T} < %.2f GeV/c;m^{2} [GeV/c^{2}]^{2};N_{counts}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -2., 5.);
+    hNsigProtonMSqrAfter[i] = new TH1F(Form("hNsigProtonMSqrAfter%i", i), Form("M^{2} for %.2f < p_{T} < %.2f GeV/c;m^{2} [GeV/c^{2}]^{2};N_{counts}", ptBin[i] + cutPtMin.at(energy), ptBin[i + 1] + cutPtMin.at(energy)), 1400, -2., 5.);
+  }
+
+  int i_part=-1;
+
+  // Loop over events
+  for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
+  {
+
+    if ((iEvent + 1) % 1000 == 0)
+    {
+      std::cout << "Working on event #[" << (iEvent + 1)
+                << "/" << events2read << "]" << std::endl;
+    }
+
+    Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
+    if (!readEvent)
+    {
+      std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
+      break;
+    }
+
+    // Retrieve femtoDst
+    StFemtoDst *dst = femtoReader->femtoDst();
+
+    // Retrieve event information
+    StFemtoEvent *event = dst->event();
+    if (!event)
+    {
+      std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
+      break;
+    }
+
+    TVector3 pVtx = event->primaryVertex();
+
+    // Event selection
+    if (!isGoodEvent(event))
+      continue;
+
+    // Track analysis
+    Int_t nTracks = dst->numberOfTracks();
+
+    // Track loop
+    for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
+    {
+
+      // Retrieve i-th femto track
+      StFemtoTrack *femtoTrack = dst->track(iTrk);
+
+      if (!femtoTrack)
+        continue;
+
+      // Must be a primary track
+      if (!femtoTrack->isPrimary())
+        continue;
+
+      //Track selection
+      if (!isGoodTrackFlow(femtoTrack, energy, pVtx))
+        continue;
+
+      // Determine pt bin for PID plots
+      int i_pt = -1;
+      for (int i = 0; i < NptBins; i++)
+      {
+        if (femtoTrack->pt() >= (ptBin[i] + cutPtMin.at(energy)) && femtoTrack->pt() < (ptBin[i + 1] + cutPtMin.at(energy)))
+        {
+          i_pt = i;
+        }
+      }
+      if (i_pt == -1)
+        continue;
+
+      // Check if track has TOF signal
+      if (femtoTrack->isTofTrack())
+      {
+        if (femtoTrack->gDCA(pVtx).Mag() >= 1.)
+          continue;
+
+        hNsigPionMSqr[i_pt]->Fill(femtoTrack->nSigmaPion(), femtoTrack->massSqr());
+        hNsigKaonMSqr[i_pt]->Fill(femtoTrack->nSigmaKaon(), femtoTrack->massSqr());
+        hNsigProtonMSqr[i_pt]->Fill(femtoTrack->nSigmaProton(), femtoTrack->massSqr());
+
+        //Loop over particle species (pi,K,p)
+        i_part = -1;
+        for (int iPID = 0; iPID < 3; iPID++)
+        {
+          //Small pT
+          //pion id
+          if (femtoTrack->pt() < 0.5 &&
+              TMath::Abs(femtoTrack->nSigmaPion()) < cutNsigPID)
+          {
+            i_part = 0;
+          }
+          // kaon id
+          if (femtoTrack->pt() < 0.5 &&
+              TMath::Abs(femtoTrack->nSigmaKaon()) < cutNsigPID)
+          {
+            i_part = 1;
+          }
+          // proton id
+          if (femtoTrack->pt() < 1. &&
+              TMath::Abs(femtoTrack->nSigmaProton()) < cutNsigPID)
+          {
+            i_part = 2;
+          }
+
+          //Large pT
+          // pion id - asymmetry cut
+          if (femtoTrack->pt() >= 0.5 &&
+              TMath::Abs(femtoTrack->nSigmaPion()) < cutNsigPID &&
+              femtoTrack->massSqr() > pionMassSqr - NpidFactor * sigMsqrPion->Eval(femtoTrack->pt()) &&
+              femtoTrack->massSqr() < pionMassSqr + NpidFactor * sigMsqrPion->Eval(femtoTrack->pt()))
+          {
+            i_part = 0;
+          }
+          // kaon id - asymmetry cut
+          if (femtoTrack->pt() >= 0.5 &&
+              TMath::Abs(femtoTrack->nSigmaKaon()) < cutNsigPID &&
+              femtoTrack->massSqr() > kaonMassSqr - NpidFactor * TMath::Abs(sigMsqrKaon->Eval(femtoTrack->pt())) &&
+              femtoTrack->massSqr() < kaonMassSqr + NpidFactor * TMath::Abs(sigMsqrKaon->Eval(femtoTrack->pt())))
+          {
+            i_part = 1;
+          }
+          // proton id
+          if (femtoTrack->pt() >= 1. &&
+              TMath::Abs(femtoTrack->nSigmaProton()) < cutNsigPID &&
+              femtoTrack->massSqr() > protMassSqr - NpidFactor * sigMsqrProton->Eval(femtoTrack->pt()) &&
+              femtoTrack->massSqr() < protMassSqr + NpidFactor * sigMsqrProton->Eval(femtoTrack->pt()))
+          {
+            i_part = 2;
+          }
+        }
+        if (i_part == -1)
+          continue;
+
+        if (i_part == 0)
+          hNsigPionMSqrAfter[i_pt]->Fill(femtoTrack->massSqr());
+        if (i_part == 1)
+          hNsigKaonMSqrAfter[i_pt]->Fill(femtoTrack->massSqr());
+        if (i_part == 2)
+          hNsigProtonMSqrAfter[i_pt]->Fill(femtoTrack->massSqr());
+
+      } //if( isTofTrack() )
+
+    } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
+
+  } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
+
+  femtoReader->Finish();
+
+  TFile *output = new TFile(outFile, "recreate");
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigPionMSqr[i]->Write();
+    hNsigKaonMSqr[i]->Write();
+    hNsigProtonMSqr[i]->Write();
+  }
+
+  for (int i = 0; i < NptBins; i++)
+  {
+    hNsigPionMSqrAfter[i]->Write();
+    hNsigKaonMSqrAfter[i]->Write();
+    hNsigProtonMSqrAfter[i]->Write();
+  }
+
+  output->Close();
+
+  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
+            << std::endl;
+  timer.Stop();
+  timer.Print();
+}
+
+Bool_t isGoodEvent(StFemtoEvent *const &event)
+{
+  if (!event)
+    return false;
+  if (event == nullptr)
+    return false;
+
+  if (event->primaryVertex().Perp() > cutVtxR)
+    return false;
+  if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
+    return false;
+
+  if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
+    return false;
+
+  return true;
+}
+
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track)
+    return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits)
+    return false;
+  if (track->nHitsPoss() <= cutNhitsPoss)
+    return false;
+  if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
+    return false;
+  if (TMath::Abs(track->eta()) >= cutEta)
+    return false;
+
+  if (track->pt() <= cutPtMin.at(_energy))
+    return false;
+  if (track->pt() > cutPtMax)
+    return false;
+  if (track->ptot() > cutPMax)
+    return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
+    return false;
+  return true;
+}
+
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track)
+    return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits)
+    return false;
+  if (track->nHitsPoss() <= cutNhitsPoss)
+    return false;
+  if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
+    return false;
+  if (TMath::Abs(track->eta()) >= cutEta)
+    return false;
+
+  if (track->pt() <= cutPtMin.at(_energy))
+    return false;
+  //if (track->pt() > cutPtMax) return false;
+  if (track->ptot() > cutPMax)
+    return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
+    return false;
+  return true;
+}
+
+Double_t GetWeight(StFemtoTrack *const &track)
+{
+  Double_t weight;
+  if (track->pt() <= cutPtWeightEP)
+  {
+    weight = track->pt();
+  }
+  else
+  {
+    weight = cutPtWeightEP;
+  }
+  return weight;
+}
+
+TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
+{
+  TVector2 qv(0., 0.);
+
+  qv.Set(TMath::Cos(_harm * track->phi()), TMath::Sin(_harm * track->phi()));
+
+  return qv;
+}
+
+Double_t AngleShift(Double_t Psi, Double_t order)
+{
+  Double_t PsiCorr = Psi;
+  if (PsiCorr > TMath::Pi() / order)
+  {
+    PsiCorr = Psi - 2. * TMath::Pi() / order;
+  }
+  if (PsiCorr < -1. * TMath::Pi() / order)
+  {
+    PsiCorr = Psi + 2. * TMath::Pi() / order;
+  }
+  return PsiCorr;
+}

+ 60 - 0
macro/mainPID.cpp

@@ -0,0 +1,60 @@
+#include <iostream>
+
+#include <TROOT.h>
+#include <TStopwatch.h>
+#include <TFile.h>
+#include <TString.h>
+
+#include "FemtoDstAnalyzer_PID.C"
+
+int main(int argc, char** argv)
+{
+  TString inFileName, outFileName;
+  TStopwatch timer;
+  if (argc < 5)
+  {
+    std::cerr << "./pid -i INPUTFILE -o OUTPUTFILE" << std::endl;
+    return 10;
+  }
+
+  for (int i=1;i<argc;i++)
+  {
+    if (std::string(argv[i]) != "-i" &&
+        std::string(argv[i]) != "-o")
+    {
+      std::cerr << "\nUnknown parameter: " << argv[i] << std::endl;
+      return 11;
+    }
+    else
+    {
+      if (std::string(argv[i]) == "-i" && i != argc-1)
+      {
+        inFileName = TString(argv[++i]);
+      }
+      if (std::string(argv[i]) == "-i" && i == argc-1)
+      {
+        std::cerr << "\nInput file was not specified!" << std::endl;
+	return 12;
+      }
+      if (std::string(argv[i]) == "-o" && i != argc-1)
+      {
+        outFileName = TString(argv[++i]);
+      }
+      if (std::string(argv[i]) == "-o" && i == argc-1)
+      {
+        std::cerr << "\nOutput file was not specified!" << std::endl;
+	return 13;
+      }
+    }
+  }
+  if (inFileName == "" || outFileName == "")
+  {
+    std::cerr << "\nInput/Output file has not been set properly!" << std::endl;
+    return 14;
+  }
+
+  FemtoDstAnalyzer_PID(inFileName.Data(),outFileName.Data());
+
+  return 0;
+}
+

+ 70 - 0
macro/mainPID_QA.cpp

@@ -0,0 +1,70 @@
+#include <iostream>
+
+#include <TROOT.h>
+#include <TStopwatch.h>
+#include <TFile.h>
+#include <TString.h>
+
+#include "FemtoDstAnalyzer_PID_QA.C"
+
+int main(int argc, char** argv)
+{
+  TString inFileName, outFileName, pidFileName;
+  TStopwatch timer;
+  if (argc < 5)
+  {
+    std::cerr << "./pidQA -i INPUTFILE -o OUTPUTFILE -pid" << std::endl;
+    return 10;
+  }
+
+  for (int i=1;i<argc;i++)
+  {
+    if (std::string(argv[i]) != "-i" &&
+        std::string(argv[i]) != "-o" &&
+	std::string(argv[i]) != "-pid")
+    {
+      std::cerr << "\nUnknown parameter: " << argv[i] << std::endl;
+      return 11;
+    }
+    else
+    {
+      if (std::string(argv[i]) == "-i" && i != argc-1)
+      {
+        inFileName = TString(argv[++i]);
+      }
+      if (std::string(argv[i]) == "-i" && i == argc-1)
+      {
+        std::cerr << "\nInput file was not specified!" << std::endl;
+	return 12;
+      }
+      if (std::string(argv[i]) == "-o" && i != argc-1)
+      {
+        outFileName = TString(argv[++i]);
+      }
+      if (std::string(argv[i]) == "-o" && i == argc-1)
+      {
+        std::cerr << "\nOutput file was not specified!" << std::endl;
+	return 13;
+      }
+      if (std::string(argv[i]) == "-pid" && i != argc-1)
+      {
+        pidFileName = TString(argv[++i]);
+      }
+      if (std::string(argv[i]) == "-pid" && i == argc-1)
+      {
+        std::cerr << "\nOutput file was not specified!" << std::endl;
+	return 13;
+      }
+    }
+  }
+  if (inFileName == "" || outFileName == "" || pidFileName == "")
+  {
+    std::cerr << "\nInput/Output file has not been set properly!" << std::endl;
+    return 14;
+  }
+
+  FemtoDstAnalyzer_PID_QA(inFileName.Data(),outFileName.Data());
+
+  return 0;
+}
+

+ 36 - 0
scripts/GenerateLists.sh

@@ -0,0 +1,36 @@
+#!/bin/bash
+
+#-Set path to the femtoDst files
+INPUT_DIR=$1
+#-Set the maximum number of files in 1 list
+NUM_DIVIDES=$2
+#-Name pattern for output lists
+NAME_PATTERN=$3
+
+#-Example usage
+# . GenerateLists.sh /mnt/pool/rhic/2/nigmatkulov/femtoDst/auau/200gev/12135/ 100
+
+CURRENT_DIR=${PWD}
+OUTPUT_DIR="../lists"
+
+mkdir -p $OUTPUT_DIR
+
+TOTAL_NUM_FILES=`ls $INPUT_DIR/*.femtoDst.root | wc -l`
+echo "Total number of DST files: ${TOTAL_NUM_FILES}"
+
+MULT=$((TOTAL_NUM_FILES / NUM_DIVIDES))
+RESIDUE=$((TOTAL_NUM_FILES % NUM_DIVIDES))
+
+echo "Mult: $MULT, Residue: $RESIDUE"
+
+for i in `seq 1 $MULT`
+do
+  OUTPUT_FILE=$OUTPUT_DIR/${NAME_PATTERN}${i}.list
+  echo $OUTPUT_FILE
+  ls $INPUT_DIR/*.femtoDst.root | head -n$((i*NUM_DIVIDES)) | tail -n$NUM_DIVIDES &> $OUTPUT_FILE
+done
+
+#-Generate last filelist for residual files
+OUTPUT_FILE=$OUTPUT_DIR/${NAME_PATTERN}$((i+1)).list
+echo $OUTPUT_FILE
+ls $INPUT_DIR/*.femtoDst.root |  tail -n$RESIDUE &> $OUTPUT_FILE

+ 12 - 0
scripts/run.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+
+INPUT_FILE=$1
+OUTPUT_FILE=$2
+LOG_FILE=$3
+
+echo "Arguments:"
+echo "Input: $INPUT_FILE"
+echo "Output: $OUTPUT_FILE"
+echo "Log: $LOG_FILE"
+
+/mnt/pool/rhic/4/parfenovpeter/STAR/Analysis/hpc_scripts/build/FemtoDstAnalyzer -i $INPUT_FILE -o $OUTPUT_FILE &> $LOG_FILE

+ 22 - 0
scripts/start.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+
+INPUT_LIST_DIR=$1
+OUTPUT_DIR=/mnt/pool/rhic/4/parfenovpeter/STAR/Analysis/hpc_scripts/OUT/`date '+%Y%m%d_%H%M%S'`
+
+QUEUE=medium
+
+mkdir -p $OUTPUT_DIR/log
+mkdir -p $OUTPUT_DIR/root
+mkdir -p $OUTPUT_DIR/sge_error
+mkdir -p $OUTPUT_DIR/sge_output
+
+ls -1 $INPUT_LIST_DIR/* | while read line
+do
+  OUTPUT_ROOT=$OUTPUT_DIR/root/$RUNID/`basename ${line%.*t}`.root
+  OUTPUT_LOG=$OUTPUT_DIR/log/$RUNID/`basename ${line%.*t}`.log
+  OUTPUT_O_SGE=$OUTPUT_DIR/sge_output/$RUNID/`basename ${line%.*t}`.out
+  OUTPUT_E_SGE=$OUTPUT_DIR/sge_error/$RUNID/`basename ${line%.*t}`.err
+
+  qsub -q $QUEUE -o $OUTPUT_O_SGE -e $OUTPUT_E_SGE run.sh -F "$line $OUTPUT_ROOT $INPUT_RECENTERING $OUTPUT_LOG"
+
+done

+ 4 - 0
set_env.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+export ST_FEMTO_DST_INC_DIR=/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent
+export ST_FEMTO_DST_BUILD_DIR=$ST_FEMTO_DST_INC_DIR