Просмотр исходного кода

Added v1, v2 as a TProfile with config parser to std::vector<float> for centrality classes (b-based)

ParfenovPeter лет назад: 2
Родитель
Сommit
461c67e202
4 измененных файлов с 287 добавлено и 114 удалено
  1. 71 26
      bin/main.cpp
  2. 2 0
      qa.cfg
  3. 196 85
      utility/Utility.cxx
  4. 18 3
      utility/Utility.h

+ 71 - 26
bin/main.cpp

@@ -92,9 +92,17 @@ int main(int argc, char **argv)
   TStopwatch timer;
   timer.Start();
 
+  Bool_t IsRead = true;
+
   if (configFileName.Length() > 0)
   {
-    qaUtility::GetInstance()->ReadConfig(configFileName);
+    IsRead = qaUtility::GetInstance()->ReadConfig(configFileName);
+  }
+
+  if (!IsRead)
+  {
+    std::cerr << "Error while reading config file. Exiting..." << std::endl;
+    return 2;
   }
 
   if (qaUtility::GetInstance()->debug)
@@ -130,6 +138,13 @@ int main(int argc, char **argv)
       std::cout << "Is_v1 = " << qaUtility::GetInstance()->Is_v1 << std::endl;
       std::cout << "Cut_v1_Event_bmin = " << qaUtility::GetInstance()->Cut_v1_Event_bmin << std::endl;
       std::cout << "Cut_v1_Event_bmax = " << qaUtility::GetInstance()->Cut_v1_Event_bmax << std::endl;
+      std::cout << "Cut_v1_Event_bCent[" << qaUtility::GetInstance()->Cut_v1_Event_bCent.size() << "] = {";
+      for (int i=0; i<(int)qaUtility::GetInstance()->Cut_v1_Event_bCent.size(); i++)
+      {
+        std::cout << qaUtility::GetInstance()->Cut_v1_Event_bCent.at(i);
+        if (i != (int)qaUtility::GetInstance()->Cut_v1_Event_bCent.size()-1) std::cout << ", ";
+        else std::cout << "};" << std::endl;
+      }
       std::cout << "Cut_v1_Particle_ptmin = " << qaUtility::GetInstance()->Cut_v1_Particle_ptmin << std::endl;
       std::cout << "Cut_v1_Particle_ptmax = " << qaUtility::GetInstance()->Cut_v1_Particle_ptmax << std::endl;
       std::cout << "Cut_v1_Particle_etamin = " << qaUtility::GetInstance()->Cut_v1_Particle_etamin << std::endl;
@@ -142,6 +157,13 @@ int main(int argc, char **argv)
       std::cout << "Is_v2 = " << qaUtility::GetInstance()->Is_v2 << std::endl;
       std::cout << "Cut_v2_Event_bmin = " << qaUtility::GetInstance()->Cut_v2_Event_bmin << std::endl;
       std::cout << "Cut_v2_Event_bmax = " << qaUtility::GetInstance()->Cut_v2_Event_bmax << std::endl;
+      std::cout << "Cut_v2_Event_bCent[" << qaUtility::GetInstance()->Cut_v2_Event_bCent.size() << "] = {";
+      for (int i=0; i<(int)qaUtility::GetInstance()->Cut_v2_Event_bCent.size(); i++)
+      {
+        std::cout << qaUtility::GetInstance()->Cut_v2_Event_bCent.at(i);
+        if (i != (int)qaUtility::GetInstance()->Cut_v2_Event_bCent.size()-1) std::cout << ", ";
+        else std::cout << "};" << std::endl;
+      }
       std::cout << "Cut_v2_Particle_ptmin = " << qaUtility::GetInstance()->Cut_v2_Particle_ptmin << std::endl;
       std::cout << "Cut_v2_Particle_ptmax = " << qaUtility::GetInstance()->Cut_v2_Particle_ptmax << std::endl;
       std::cout << "Cut_v2_Particle_etamin = " << qaUtility::GetInstance()->Cut_v2_Particle_etamin << std::endl;
@@ -215,11 +237,11 @@ int main(int argc, char **argv)
   TH1F *h_refmult_Particle_PID_y[qaUtility::GetInstance()->npid];
   TH1F *h_refmult_Particle_PID_z[qaUtility::GetInstance()->npid];
 
-  TProfile2D *p_v1_PID_pt_b[qaUtility::GetInstance()->npid];
-  TProfile2D *p_v1_PID_y_b[qaUtility::GetInstance()->npid];
+  TProfile *p_v1_PID_pt[qaUtility::GetInstance()->maxCentBins][qaUtility::GetInstance()->npid];
+  TProfile *p_v1_PID_y[qaUtility::GetInstance()->maxCentBins][qaUtility::GetInstance()->npid];
 
-  TProfile2D *p_v2_PID_pt_b[qaUtility::GetInstance()->npid];
-  TProfile2D *p_v2_PID_y_b[qaUtility::GetInstance()->npid];
+  TProfile *p_v2_PID_pt[qaUtility::GetInstance()->maxCentBins][qaUtility::GetInstance()->npid];
+  TProfile *p_v2_PID_y[qaUtility::GetInstance()->maxCentBins][qaUtility::GetInstance()->npid];
 
   for (int i=0; i<qaUtility::GetInstance()->npid; i++)
   {
@@ -246,12 +268,17 @@ int main(int argc, char **argv)
     h_refmult_Particle_PID_x[i] = new TH1F(Form("h_refmult_Particle_PID_x_%i",i),Form("h_refmult_Particle_PID_x_%i;x, fm; dN/dx",i), 800, -8., 8.);
     h_refmult_Particle_PID_y[i] = new TH1F(Form("h_refmult_Particle_PID_y_%i",i),Form("h_refmult_Particle_PID_y_%i;y, fm; dN/dy",i), 800, -8., 8.);
     h_refmult_Particle_PID_z[i] = new TH1F(Form("h_refmult_Particle_PID_z_%i",i),Form("h_refmult_Particle_PID_z_%i;z, fm; dN/dz",i), 800, -8., 8.);
+  }
+  for (int i=0; i<qaUtility::GetInstance()->maxCentBins; i++)
+  {
+    for (int j=0; j<qaUtility::GetInstance()->npid; j++)
+    {
+      p_v1_PID_pt[i][j] = new TProfile(Form("p_v1_PID_pt_%i_%i",i,j),Form("p_v1_PID_pt_%i_%i;p_{T}, GeV/c;v_{1}",i,j),100,0.,5.);
+      p_v1_PID_y[i][j] = new TProfile(Form("p_v1_PID_y_%i_%i",i,j),Form("p_v1_PID_y_%i_%i;y;v_{1}",i,j),400,-10.,10.);
 
-    p_v1_PID_pt_b[i] = new TProfile2D(Form("p_v1_PID_pt_b_%i",i),Form("p_v1_PID_pt_b_%i;p_{T}, GeV/c;b, fm;v_{1}",i),500,0.,5.,2000,0.,20.);
-    p_v1_PID_y_b[i] = new TProfile2D(Form("p_v1_PID_y_b_%i",i),Form("p_v1_PID_y_b_%i;y;b, fm;v_{1}",i),400,-10.,10.,2000,0.,20.);
-
-    p_v2_PID_pt_b[i] = new TProfile2D(Form("p_v2_PID_pt_b_%i",i),Form("p_v2_PID_pt_b_%i;p_{T}, GeV/c;b, fm;v_{1}",i),500,0.,5.,2000,0.,20.);
-    p_v2_PID_y_b[i] = new TProfile2D(Form("p_v2_PID_y_b_%i",i),Form("p_v2_PID_y_b_%i;y;b, fm;v_{1}",i),400,-10.,10.,2000,0.,20.);
+      p_v2_PID_pt[i][j] = new TProfile(Form("p_v2_PID_pt_%i_%i",i,j),Form("p_v2_PID_pt_%i_%i;p_{T}, GeV/c;v_{2}",i,j),100,0.,5.);
+      p_v2_PID_y[i][j] = new TProfile(Form("p_v2_PID_y_%i_%i",i,j),Form("p_v2_PID_y_%i_%i;y;v_{2}",i,j),400,-10.,10.);
+    }
   }
 
   qaReader_manager *readerManager;
@@ -293,7 +320,7 @@ int main(int argc, char **argv)
   Long64_t Minbias_counter = 0;
   Int_t ipid;
   Double_t v1, v2, y;
-  Int_t eta_w;
+  Int_t eta_w, centBin;
   while(Minbias_counter < Nentries)
   {
     if (Minbias_counter % 1000 == 0)
@@ -434,12 +461,12 @@ int main(int argc, char **argv)
       }
 
       if (qaUtility::GetInstance()->Cut_Event_v1(event) && 
-        qaUtility::GetInstance()->Cut_Particle_v1(particle) &&
         qaUtility::GetInstance()->Is_v1 == 1)
       {
         v1 = TMath::Cos(1.*(particle->GetPhi() - event->GetPhiRP()));
-        eta_w = (particle->GetEta() >= 0.) ? 1 : -1;
+        eta_w = 1; //(particle->GetEta() >= 0.) ? 1 : -1;
         y = 0.5*TMath::Log( (particle->GetEnergy() + particle->GetPz())/(particle->GetEnergy() - particle->GetPz()) );
+        centBin = qaUtility::GetInstance()->GetCentralityBin(event->GetB(), qaUtility::GetInstance()->Cut_v1_Event_bCent);
 
         if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) > 0)
         {
@@ -450,19 +477,25 @@ int main(int argc, char **argv)
           ipid = 4;
         }
         ipid = qaUtility::GetInstance()->GetPdgId(particle->GetPdg());
-        if (ipid != -1 && ipid != 0 && ipid != 4)
+        if (ipid != -1 && ipid != 0 && ipid != 4 && centBin != -1)
         {
-          p_v1_PID_pt_b[ipid]->Fill(particle->GetPt(),event->GetB(),v1*eta_w);
-          p_v1_PID_y_b[ipid]->Fill(y,event->GetB(),v1);
+          if (qaUtility::GetInstance()->Cut_Particle_v1_pt(particle))
+          {
+            p_v1_PID_pt[centBin][ipid]->Fill(particle->GetPt(),v1*eta_w);
+          }
+          if (qaUtility::GetInstance()->Cut_Particle_v1_y(particle))
+          {
+            p_v1_PID_y[centBin][ipid]->Fill(y,v1);
+          }
         }
       }
 
       if (qaUtility::GetInstance()->Cut_Event_v2(event) && 
-        qaUtility::GetInstance()->Cut_Particle_v2(particle) &&
         qaUtility::GetInstance()->Is_v2 == 1)
       {
         v2 = TMath::Cos(2.*(particle->GetPhi() - event->GetPhiRP()));
         y = 0.5*TMath::Log( (particle->GetEnergy() + particle->GetPz())/(particle->GetEnergy() - particle->GetPz()) );
+        centBin = qaUtility::GetInstance()->GetCentralityBin(event->GetB(), qaUtility::GetInstance()->Cut_v2_Event_bCent);
 
         if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) > 0)
         {
@@ -473,10 +506,16 @@ int main(int argc, char **argv)
           ipid = 4;
         }
         ipid = qaUtility::GetInstance()->GetPdgId(particle->GetPdg());
-        if (ipid != -1 && ipid != 0 && ipid != 4)
+        if (ipid != -1 && ipid != 0 && ipid != 4 && centBin != -1)
         {
-          p_v2_PID_pt_b[ipid]->Fill(particle->GetPt(),event->GetB(),v2);
-          p_v2_PID_y_b[ipid]->Fill(y,event->GetB(),v2);
+          if (qaUtility::GetInstance()->Cut_Particle_v2_pt(particle))
+          {
+            p_v2_PID_pt[centBin][ipid]->Fill(particle->GetPt(),v2);
+          }
+          if (qaUtility::GetInstance()->Cut_Particle_v2_y(particle))
+          {
+            p_v2_PID_y[centBin][ipid]->Fill(y,v2);
+          }
         }
       }
 
@@ -597,10 +636,13 @@ int main(int argc, char **argv)
     fo->mkdir("v1");
     fo->cd("v1");
 
-    for (int i=0; i<qaUtility::GetInstance()->npid; i++)
+    for (int i=0; i<(int)qaUtility::GetInstance()->Cut_v1_Event_bCent.size()-1; i++)
     {
-      p_v1_PID_pt_b[i]->Write();
-      p_v1_PID_y_b[i]->Write();
+      for (int j=0; j<qaUtility::GetInstance()->npid; j++)
+      {
+        p_v1_PID_pt[i][j]->Write();
+        p_v1_PID_y[i][j]->Write();
+      }
     }
   }
 
@@ -609,10 +651,13 @@ int main(int argc, char **argv)
     fo->mkdir("v2");
     fo->cd("v2");
 
-    for (int i=0; i<qaUtility::GetInstance()->npid; i++)
+    for (int i=0; i<(int)qaUtility::GetInstance()->Cut_v2_Event_bCent.size()-1; i++)
     {
-      p_v2_PID_pt_b[i]->Write();
-      p_v2_PID_y_b[i]->Write();
+      for (int j=0; j<qaUtility::GetInstance()->npid; j++)
+      {
+        p_v2_PID_pt[i][j]->Write();
+        p_v2_PID_y[i][j]->Write();
+      }
     }
   }
 

+ 2 - 0
qa.cfg

@@ -19,6 +19,7 @@ Cut_refmult_Particle_etamax:          0.5
 Is_v1                                 1
 Cut_v1_Event_bmin:                    0.
 Cut_v1_Event_bmax:                    16.
+Cut_v1_Event_bCent:                   6.6 8.1
 Cut_v1_Particle_ptmin:                1.0
 Cut_v1_Particle_ptmax:                1.5
 Cut_v1_Particle_etamin:               -100.
@@ -28,6 +29,7 @@ Cut_v1_Particle_ymax:                 -0.15
 Is_v2                                 1
 Cut_v2_Event_bmin:                    0.
 Cut_v2_Event_bmax:                    16.
+Cut_v2_Event_bCent:                   6.6 8.1
 Cut_v2_Particle_ptmin:                1.0
 Cut_v2_Particle_ptmax:                1.5
 Cut_v2_Particle_etamin:               -100.

+ 196 - 85
utility/Utility.cxx

@@ -4,45 +4,50 @@
 
 ClassImp(qaUtility);
 
-qaUtility* qaUtility::fUtility = nullptr;;
+qaUtility *qaUtility::fUtility = nullptr;
+;
 
 qaUtility::qaUtility() : Nevents(-1),
-debug(0),
-format("mctree"),
-Is_minbias(1),
-Is_refmult(1),
-Is_v1(1),
-Is_v2(1),
-Cut_minbias_Event_bmin(0.),
-Cut_minbias_Event_bmax(16.),
-Cut_minbias_Particle_ptmin(0.),
-Cut_minbias_Particle_ptmax(100.),
-Cut_minbias_Particle_etamin(-100.),
-Cut_minbias_Particle_etamax(100.),
-Cut_minbias_Particle_ymin(-100.),
-Cut_minbias_Particle_ymax(100.),
-Cut_refmult_Event_bmin(0.),
-Cut_refmult_Event_bmax(16.),
-Cut_refmult_Particle_ptmin(0.15),
-Cut_refmult_Particle_ptmax(3.),
-Cut_refmult_Particle_etamin(-0.5),
-Cut_refmult_Particle_etamax(0.5),
-Cut_v1_Event_bmin(0.),
-Cut_v1_Event_bmax(16.),
-Cut_v1_Particle_ptmin(1.0),
-Cut_v1_Particle_ptmax(1.5),
-Cut_v1_Particle_etamin(-100.),
-Cut_v1_Particle_etamax(100.),
-Cut_v1_Particle_ymin(-0.25),
-Cut_v1_Particle_ymax(-0.15),
-Cut_v2_Event_bmin(0.),
-Cut_v2_Event_bmax(16.),
-Cut_v2_Particle_ptmin(1.0),
-Cut_v2_Particle_ptmax(1.5),
-Cut_v2_Particle_etamin(-100.),
-Cut_v2_Particle_etamax(100.),
-Cut_v2_Particle_ymin(-0.05),
-Cut_v2_Particle_ymax(0.05)
+                         debug(0),
+                         format("mctree"),
+                         Is_minbias(1),
+                         Is_refmult(1),
+                         Is_v1(1),
+                         Is_v2(1),
+                         Cut_minbias_Event_bmin(0.),
+                         Cut_minbias_Event_bmax(16.),
+                         Cut_minbias_Particle_ptmin(0.),
+                         Cut_minbias_Particle_ptmax(100.),
+                         Cut_minbias_Particle_etamin(-100.),
+                         Cut_minbias_Particle_etamax(100.),
+                         Cut_minbias_Particle_ymin(-100.),
+                         Cut_minbias_Particle_ymax(100.),
+                         Cut_refmult_Event_bmin(0.),
+                         Cut_refmult_Event_bmax(16.),
+                         Cut_refmult_Particle_ptmin(0.15),
+                         Cut_refmult_Particle_ptmax(3.),
+                         Cut_refmult_Particle_etamin(-0.5),
+                         Cut_refmult_Particle_etamax(0.5),
+                         Cut_v1_Event_bmin(0.),
+                         Cut_v1_Event_bmax(16.),
+                         sCut_v1_Event_bCent(""),
+                         Cut_v1_Event_bCent({}),
+                         Cut_v1_Particle_ptmin(1.0),
+                         Cut_v1_Particle_ptmax(1.5),
+                         Cut_v1_Particle_etamin(-100.),
+                         Cut_v1_Particle_etamax(100.),
+                         Cut_v1_Particle_ymin(-0.25),
+                         Cut_v1_Particle_ymax(-0.15),
+                         Cut_v2_Event_bmin(0.),
+                         Cut_v2_Event_bmax(16.),
+                         sCut_v2_Event_bCent(""),
+                         Cut_v2_Event_bCent({}),
+                         Cut_v2_Particle_ptmin(1.0),
+                         Cut_v2_Particle_ptmax(1.5),
+                         Cut_v2_Particle_etamin(-100.),
+                         Cut_v2_Particle_etamax(100.),
+                         Cut_v2_Particle_ymin(-0.05),
+                         Cut_v2_Particle_ymax(0.05)
 {
 }
 
@@ -94,6 +99,7 @@ Bool_t qaUtility::ReadConfig(const TString &configFileName)
   Is_v1 = env.GetValue("Is_v1", 0);
   Cut_v1_Event_bmin = env.GetValue("Cut_v1_Event_bmin", 0.);
   Cut_v1_Event_bmax = env.GetValue("Cut_v1_Event_bmax", 0.);
+  sCut_v1_Event_bCent = env.GetValue("Cut_v1_Event_bCent", "");
 
   Cut_v1_Particle_ptmin = env.GetValue("Cut_v1_Particle_ptmin", 0.);
   Cut_v1_Particle_ptmax = env.GetValue("Cut_v1_Particle_ptmax", 0.);
@@ -105,6 +111,7 @@ Bool_t qaUtility::ReadConfig(const TString &configFileName)
   Is_v2 = env.GetValue("Is_v2", 0);
   Cut_v2_Event_bmin = env.GetValue("Cut_v2_Event_bmin", 0.);
   Cut_v2_Event_bmax = env.GetValue("Cut_v2_Event_bmax", 0.);
+  sCut_v2_Event_bCent = env.GetValue("Cut_v2_Event_bCent", "");
 
   Cut_v2_Particle_ptmin = env.GetValue("Cut_v2_Particle_ptmin", 0.);
   Cut_v2_Particle_ptmax = env.GetValue("Cut_v2_Particle_ptmax", 0.);
@@ -113,6 +120,8 @@ Bool_t qaUtility::ReadConfig(const TString &configFileName)
   Cut_v2_Particle_ymin = env.GetValue("Cut_v2_Particle_ymin", 0.);
   Cut_v2_Particle_ymax = env.GetValue("Cut_v2_Particle_ymax", 0.);
 
+  if (!initCentrality()) return false;
+
   return true;
 }
 
@@ -129,104 +138,206 @@ TChain *qaUtility::initChain(const TString &inputFileName, const char *chainName
   return chain;
 }
 
-Bool_t qaUtility::Cut_Event_minbias(qaEvent *const& event)
+std::vector<Float_t> qaUtility::ReadBvector(std::string _input)
+{
+  std::vector<Float_t> vB;
+
+  std::istringstream iss(_input);
+
+  std::copy(std::istream_iterator<Float_t>(iss),
+            std::istream_iterator<Float_t>(),
+            std::back_inserter(vB));
+
+  return vB;
+}
+
+Bool_t qaUtility::initCentrality()
 {
-  if (event->GetB() < Cut_minbias_Event_bmin) return false;
-  if (event->GetB() > Cut_minbias_Event_bmax) return false;
+  Cut_v1_Event_bCent = ReadBvector(sCut_v1_Event_bCent);
+  Cut_v2_Event_bCent = ReadBvector(sCut_v2_Event_bCent);
 
+  if (Cut_v1_Event_bCent.size() == 0 || 
+    Cut_v2_Event_bCent.size() == 0)
+  {
+    return false;
+  }
+  
   return true;
 }
 
-Bool_t qaUtility::Cut_Event_refmult(qaEvent *const& event)
+Bool_t qaUtility::Cut_Event_minbias(qaEvent *const &event)
 {
-  if (event->GetB() < Cut_refmult_Event_bmin) return false;
-  if (event->GetB() > Cut_refmult_Event_bmax) return false;
+  if (event->GetB() < Cut_minbias_Event_bmin)
+    return false;
+  if (event->GetB() > Cut_minbias_Event_bmax)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Event_v1(qaEvent *const& event)
+Bool_t qaUtility::Cut_Event_refmult(qaEvent *const &event)
 {
-  if (event->GetB() < Cut_v1_Event_bmin) return false;
-  if (event->GetB() > Cut_v1_Event_bmax) return false;
+  if (event->GetB() < Cut_refmult_Event_bmin)
+    return false;
+  if (event->GetB() > Cut_refmult_Event_bmax)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Event_v2(qaEvent *const& event)
+Bool_t qaUtility::Cut_Event_v1(qaEvent *const &event)
 {
-  if (event->GetB() < Cut_v2_Event_bmin) return false;
-  if (event->GetB() > Cut_v2_Event_bmax) return false;
+  if (event->GetB() < Cut_v1_Event_bmin)
+    return false;
+  if (event->GetB() > Cut_v1_Event_bmax)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Particle_minbias(qaParticle *const& particle)
+Bool_t qaUtility::Cut_Event_v2(qaEvent *const &event)
 {
-  if (particle->GetPt() < Cut_minbias_Particle_ptmin) return false;
-  if (particle->GetPt() > Cut_minbias_Particle_ptmax) return false;
-  if (particle->GetEta() < Cut_minbias_Particle_etamin) return false;
-  if (particle->GetEta() > Cut_minbias_Particle_etamax) return false;
+  if (event->GetB() < Cut_v2_Event_bmin)
+    return false;
+  if (event->GetB() > Cut_v2_Event_bmax)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Particle_refmult(qaParticle *const& particle)
+Bool_t qaUtility::Cut_Particle_minbias(qaParticle *const &particle)
 {
-  if (particle->GetPt() < Cut_refmult_Particle_ptmin) return false;
-  if (particle->GetPt() > Cut_refmult_Particle_ptmax) return false;
-  if (particle->GetEta() < Cut_refmult_Particle_etamin) return false;
-  if (particle->GetEta() > Cut_refmult_Particle_etamax) return false;
-  
+  if (particle->GetPt() < Cut_minbias_Particle_ptmin)
+    return false;
+  if (particle->GetPt() > Cut_minbias_Particle_ptmax)
+    return false;
+  if (particle->GetEta() < Cut_minbias_Particle_etamin)
+    return false;
+  if (particle->GetEta() > Cut_minbias_Particle_etamax)
+    return false;
+
+  return true;
+}
+
+Bool_t qaUtility::Cut_Particle_refmult(qaParticle *const &particle)
+{
+  if (particle->GetPt() < Cut_refmult_Particle_ptmin)
+    return false;
+  if (particle->GetPt() > Cut_refmult_Particle_ptmax)
+    return false;
+  if (particle->GetEta() < Cut_refmult_Particle_etamin)
+    return false;
+  if (particle->GetEta() > Cut_refmult_Particle_etamax)
+    return false;
+
   Double_t charge = GetCharge(particle->GetPdg());
-  if (charge == error_code) return false;
-  if (charge == 0) return false;
+  if (charge == error_code)
+    return false;
+  if (charge == 0)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Particle_v1(qaParticle *const& particle)
+Bool_t qaUtility::Cut_Particle_v1_pt(qaParticle *const &particle)
 {
-  if (particle->GetPt() < Cut_v1_Particle_ptmin) return false;
-  if (particle->GetPt() > Cut_v1_Particle_ptmax) return false;
-  if (particle->GetEta() < Cut_v1_Particle_etamin) return false;
-  if (particle->GetEta() > Cut_v1_Particle_etamax) return false;
+  double y = 0.5 * TMath::Log( (particle->GetEnergy() + particle->GetPz())/(particle->GetEnergy() - particle->GetPz()) );
+  if (y < Cut_v1_Particle_ymin)
+    return false;
+  if (y > Cut_v1_Particle_ymax)
+    return false;
+
+  Double_t charge = GetCharge(particle->GetPdg());
+  if (charge == error_code)
+    return false;
+  if (charge == 0)
+    return false;
+
+  return true;
+}
+
+Bool_t qaUtility::Cut_Particle_v1_y(qaParticle *const &particle)
+{
+  if (particle->GetPt() < Cut_v1_Particle_ptmin)
+    return false;
+  if (particle->GetPt() > Cut_v1_Particle_ptmax)
+    return false;
+
+  Double_t charge = GetCharge(particle->GetPdg());
+  if (charge == error_code)
+    return false;
+  if (charge == 0)
+    return false;
+
+  return true;
+}
+
+Bool_t qaUtility::Cut_Particle_v2_pt(qaParticle *const &particle)
+{
+  double y = 0.5 * TMath::Log( (particle->GetEnergy() + particle->GetPz())/(particle->GetEnergy() - particle->GetPz()) );
+  if (y < Cut_v2_Particle_ymin)
+    return false;
+  if (y > Cut_v2_Particle_ymax)
+    return false;
 
   Double_t charge = GetCharge(particle->GetPdg());
-  if (charge == error_code) return false;
-  if (charge == 0) return false;
+  if (charge == error_code)
+    return false;
+  if (charge == 0)
+    return false;
 
   return true;
 }
 
-Bool_t qaUtility::Cut_Particle_v2(qaParticle *const& particle)
+Bool_t qaUtility::Cut_Particle_v2_y(qaParticle *const &particle)
 {
-  if (particle->GetPt() < Cut_v2_Particle_ptmin) return false;
-  if (particle->GetPt() > Cut_v2_Particle_ptmax) return false;
-  if (particle->GetEta() < Cut_v2_Particle_etamin) return false;
-  if (particle->GetEta() > Cut_v2_Particle_etamax) return false;
+  if (particle->GetPt() < Cut_v2_Particle_ptmin)
+    return false;
+  if (particle->GetPt() > Cut_v2_Particle_ptmax)
+    return false;
 
   Double_t charge = GetCharge(particle->GetPdg());
-  if (charge == error_code) return false;
-  if (charge == 0) return false;
+  if (charge == error_code)
+    return false;
+  if (charge == 0)
+    return false;
 
   return true;
 }
 
 Double_t qaUtility::GetCharge(Int_t pdg)
 {
-  auto particle = (TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg);
-  if (!particle) return error_code;
-  return particle->Charge()/3.;
+  auto particle = (TParticlePDG *)TDatabasePDG::Instance()->GetParticle(pdg);
+  if (!particle)
+    return error_code;
+  return particle->Charge() / 3.;
 }
 
 Int_t qaUtility::GetPdgId(Int_t pdg)
 {
-  if (pdg == vpdg.at(1)) return 1;
-  if (pdg == vpdg.at(2)) return 2;
-  if (pdg == vpdg.at(3)) return 3;
-  if (pdg == vpdg.at(5)) return 5;
-  if (pdg == vpdg.at(6)) return 6;
-  if (pdg == vpdg.at(7)) return 7;
+  if (pdg == vpdg.at(1))
+    return 1;
+  if (pdg == vpdg.at(2))
+    return 2;
+  if (pdg == vpdg.at(3))
+    return 3;
+  if (pdg == vpdg.at(5))
+    return 5;
+  if (pdg == vpdg.at(6))
+    return 6;
+  if (pdg == vpdg.at(7))
+    return 7;
   return -1;
 }
+
+Int_t qaUtility::GetCentralityBin(Float_t b, std::vector<Float_t> vcent)
+{
+  if (vcent.size() == 0) return -1;
+
+  for (int i=0; i<(int)vcent.size()-1; i++)
+  {
+    if (b >= vcent.at(i) && b < vcent.at(i+1)) return i;
+  }
+
+  return -1;
+}

+ 18 - 3
utility/Utility.h

@@ -4,6 +4,11 @@
 #include <iostream>
 #include <fstream>
 #include <vector>
+#include <algorithm>
+#include <iterator>
+#include <cassert>
+#include <sstream>
+#include <string>
 
 #include <Rtypes.h>
 #include <TString.h>
@@ -30,6 +35,7 @@ public:
   const Int_t npid = 8;
   const std::vector<Int_t> vpdg = {0, 211, 321, 2212, 0, -211, -321, -2212};
   const std::vector<Double_t> mpdg = {error_code, 0.13957, 0.493677, 0.938272, error_code, 0.13957, 0.493677, 0.938272};
+  const Int_t maxCentBins = 20;
 
   Int_t Nevents;
   Int_t debug;
@@ -58,6 +64,8 @@ public:
 
   Double_t Cut_v1_Event_bmin;
   Double_t Cut_v1_Event_bmax;
+  std::string sCut_v1_Event_bCent;
+  std::vector<Float_t> Cut_v1_Event_bCent;
   Double_t Cut_v1_Particle_ptmin;
   Double_t Cut_v1_Particle_ptmax;
   Double_t Cut_v1_Particle_etamin;
@@ -67,6 +75,8 @@ public:
 
   Double_t Cut_v2_Event_bmin;
   Double_t Cut_v2_Event_bmax;
+  std::string sCut_v2_Event_bCent;
+  std::vector<Float_t> Cut_v2_Event_bCent;
   Double_t Cut_v2_Particle_ptmin;
   Double_t Cut_v2_Particle_ptmax;
   Double_t Cut_v2_Particle_etamin;
@@ -77,19 +87,24 @@ public:
   Bool_t ReadConfig(const TString &configFileName);
   TChain *initChain(const TString &inputFileName, const char *chainName);
 
+  std::vector<Float_t> ReadBvector(std::string _input);
+  Bool_t initCentrality();
+
   Bool_t Cut_Event_minbias(qaEvent *const &event);
   Bool_t Cut_Event_refmult(qaEvent *const &event);
   Bool_t Cut_Event_v1(qaEvent *const &event);
   Bool_t Cut_Event_v2(qaEvent *const &event);
   Bool_t Cut_Particle_minbias(qaParticle *const &particle);
   Bool_t Cut_Particle_refmult(qaParticle *const &particle);
-  Bool_t Cut_Particle_v1(qaParticle *const &particle);
-  Bool_t Cut_Particle_v2(qaParticle *const &particle);
+  Bool_t Cut_Particle_v1_pt(qaParticle *const &particle);
+  Bool_t Cut_Particle_v1_y(qaParticle *const &particle);
+  Bool_t Cut_Particle_v2_pt(qaParticle *const &particle);
+  Bool_t Cut_Particle_v2_y(qaParticle *const &particle);
 
   Double_t GetCharge(Int_t pdg);
   Int_t GetPdgId(Int_t pdg);
 
-  
+  Int_t GetCentralityBin(Float_t b, std::vector<Float_t> vcent);
 
   ClassDef(qaUtility,0);
 }; // class qaUtility