Reworked converter

CMakeLists.txt

@@ -47,73 +47,70 @@ set(CMAKE_BUILD_TYPE Debug)
 # set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -ffast-math")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall")
 
-set(MCONVERTER_INCLUDE_DIRECTORIES
+set(QATOOLS_INCLUDE_DIRECTORIES
   ${CMAKE_CURRENT_SOURCE_DIR}
   ${CMAKE_CURRENT_SOURCE_DIR}/bin
   ${CMAKE_CURRENT_SOURCE_DIR}/format
   ${CMAKE_CURRENT_SOURCE_DIR}/readers
-  ${CMAKE_CURRENT_SOURCE_DIR}/writers
   ${CMAKE_CURRENT_SOURCE_DIR}/utility
   ${ROOT_INLCUDE_DIRS}
 )
 
-set(MCONVERTER_INCLUDE_LIBRARIES
+set(QATOOLS_INCLUDE_LIBRARIES
   ${ROOT_LIBRARIES}
 )
 
-set(MCONVERTER_LIBRARY_h_files
-  ${CMAKE_CURRENT_SOURCE_DIR}/format/mciParticle.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/format/mciEvent.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_smash_root.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_mcpico.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_manager.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/writers/mciWriter_mcpico.h
-  ${CMAKE_CURRENT_SOURCE_DIR}/writers/mciWriter_manager.h
+set(QATOOLS_LIBRARY_h_files
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaParticle.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaParticleLight.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaEvent.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_smash_root.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_mcpico.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_manager.h
   ${CMAKE_CURRENT_SOURCE_DIR}/utility/Utility.h
 )
 
-set(MCONVERTER_LIBRARY_cxx_files
-  ${CMAKE_CURRENT_SOURCE_DIR}/format/mciParticle.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/format/mciEvent.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_smash_root.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_mcpico.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_manager.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/writers/mciWriter_mcpico.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/writers/mciWriter_manager.cxx
+set(QATOOLS_LIBRARY_cxx_files
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaParticle.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaParticleLight.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/format/qaEvent.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_smash_root.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_mcpico.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_manager.cxx
   ${CMAKE_CURRENT_SOURCE_DIR}/utility/Utility.cxx
 )
 
 if(MCINI)
-  list(APPEND MCONVERTER_INCLUDE_DIRECTORIES $ENV{MCINI}/include)
-  list(APPEND MCONVERTER_INCLUDE_LIBRARIES ${MCINI})
-  list(APPEND MCONVERTER_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_mcini.h)
-  list(APPEND MCONVERTER_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_mcini.cxx)
+  list(APPEND QATOOLS_INCLUDE_DIRECTORIES $ENV{MCINI}/include)
+  list(APPEND QATOOLS_INCLUDE_LIBRARIES ${MCINI})
+  list(APPEND QATOOLS_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_mcini.h)
+  list(APPEND QATOOLS_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_mcini.cxx)
 endif()
 
 if(PHQMD)
-  list(APPEND MCONVERTER_INCLUDE_DIRECTORIES $ENV{PHQMD_PATH})
-  list(APPEND MCONVERTER_INCLUDE_LIBRARIES ${PHQMD})
-  list(APPEND MCONVERTER_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_phqmd.h)
-  list(APPEND MCONVERTER_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_phqmd.cxx)
-  list(APPEND MCONVERTER_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_hsd_root.h)
-  list(APPEND MCONVERTER_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/mciReader_hsd_root.cxx)
+  list(APPEND QATOOLS_INCLUDE_DIRECTORIES $ENV{PHQMD_PATH})
+  list(APPEND QATOOLS_INCLUDE_LIBRARIES ${PHQMD})
+  list(APPEND QATOOLS_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_phqmd.h)
+  list(APPEND QATOOLS_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_phqmd.cxx)
+  list(APPEND QATOOLS_LIBRARY_h_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_hsd_root.h)
+  list(APPEND QATOOLS_LIBRARY_cxx_files ${CMAKE_CURRENT_SOURCE_DIR}/readers/qaReader_hsd_root.cxx)
 endif()
 
-include_directories(${MCONVERTER_INCLUDE_DIRECTORIES})
+include_directories(${QATOOLS_INCLUDE_DIRECTORIES})
 
-set(MCONVERTER_LinkDef
-  ${CMAKE_CURRENT_SOURCE_DIR}/mConverter.LinkDef.h
+set(QATOOLS_LinkDef
+  ${CMAKE_CURRENT_SOURCE_DIR}/qaTools.LinkDef.h
 )
 
 #---Generate dictionary
-ROOT_GENERATE_DICTIONARY(G__mConverter
-  ${MCONVERTER_LIBRARY_h_files}
-  LINKDEF ${MCONVERTER_LinkDef}
+ROOT_GENERATE_DICTIONARY(G__qaTools
+  ${QATOOLS_LIBRARY_h_files}
+  LINKDEF ${QATOOLS_LinkDef}
 )
 
 #---Compile library
-add_library(mConverter SHARED ${MCONVERTER_LIBRARY_cxx_files} G__mConverter.cxx)
-target_link_libraries(mConverter ${MCONVERTER_INCLUDE_LIBRARIES})
+add_library(qaTools SHARED ${QATOOLS_LIBRARY_cxx_files} G__qaTools.cxx)
+target_link_libraries(qaTools ${QATOOLS_INCLUDE_LIBRARIES})
 
 # Get the current working branch
 execute_process(
@@ -151,11 +148,11 @@ if( GIT_TAG_VERSION_ERROR_CODE )
   set(GIT_TAG_VERSION 0.0)
 endif()
 
-set_target_properties (mConverter 
+set_target_properties (qaTools 
   PROPERTIES 
   VERSION ${GIT_TAG_VERSION} SOVERSION ${GIT_BRANCH}-${GIT_COMMIT_HASH}
 )
 
 #---Compile main executable
-add_executable(ModelConverter "${CMAKE_CURRENT_SOURCE_DIR}/bin/main.cpp")
-target_link_libraries(ModelConverter mConverter)
+add_executable(mConvert "${CMAKE_CURRENT_SOURCE_DIR}/bin/main.cpp")
+target_link_libraries(mConvert qaTools)

bin/main.cpp

@@ -2,36 +2,40 @@
 #include <vector>
 
 #include <TString.h>
+#include <TH1D.h>
+#include <TH1I.h>
+#include <TH2D.h>
+#include <TProfile.h>
+#include <TProfile3D.h>
 #include <TFile.h>
 #include <TStopwatch.h>
 
-#include <mciParticle.h>
-#include <mciEvent.h>
-#include <mciReader_manager.h>
-#include <mciReader_smash_root.h>
-#include <mciReader_mcpico.h>
-#include <mciWriter_manager.h>
-#include <mciWriter_mcpico.h>
+#include <qaParticle.h>
+#include <qaParticleLight.h>
+#include <qaEvent.h>
+#include <qaReader_manager.h>
+#include <qaReader_smash_root.h>
+#include <qaReader_mcpico.h>
 #include <Utility.h>
 
 #ifdef _MCINI_
-#include <mciReader_mcini.h>
+#include <qaReader_mcini.h>
 #endif
 #ifdef _PHQMD_
-#include <mciReader_phqmd.h>
+#include <qaReader_phqmd.h>
 #endif
 #ifdef _HSD_ROOT_
-#include <mciReader_hsd_root.h>
+#include <qaReader_hsd_root.h>
 #endif
 
 int main(int argc, char **argv)
 {
-  TString iFileName, oFileName;
-  
-  if (argc < 9)
+  TString iFileName, oFileName, configFileName = "";
+
+  if (argc < 7)
   {
-    std::cerr << "./ModelConverter -i input.list -o output.root -input-format [FORMAT] -output-format [FORMAT] [OPTIONAL: -debug]" << std::endl;
-    std::cerr << "Available input formats:" << std::endl;
+    std::cerr << "./mConvert -i input.list -o qa_output.root -format [FORMAT] [OPTIONAL: -config qa.cfg]" << std::endl;
+    std::cerr << "Available formats:" << std::endl;
     std::cerr << "\tmcpico   - simple custom ROOT format to store model data." << std::endl;
     std::cerr << "\tparticle - ROOT format that is used by the SMASH model." << std::endl;
   #ifdef _MCINI_
@@ -43,17 +47,14 @@ int main(int argc, char **argv)
   #ifdef _HSD_ROOT_
     std::cerr << "\thsd      - custom ROOT format to store HSD model data." << std::endl;
   #endif
-  std::cerr << "Available output formats:" << std::endl;
-    std::cerr << "\tmcpico   - simple custom ROOT format to store model data." << std::endl;
     return 1;
   }
   for (int i = 1; i < argc; i++)
   {
     if (std::string(argv[i]) != "-i" &&
         std::string(argv[i]) != "-o" &&
-        std::string(argv[i]) != "-input-format" &&
-        std::string(argv[i]) != "-output-format" &&
-        std::string(argv[i]) != "-debug")
+        std::string(argv[i]) != "-format" &&
+        std::string(argv[i]) != "-config")
     {
       std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
       return 2;
@@ -80,72 +81,106 @@ int main(int argc, char **argv)
         std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
         return 1;
       }
-      if (std::string(argv[i]) == "-input-format" && i != argc - 1)
+      if (std::string(argv[i]) == "-format" && i != argc - 1)
       {
-        mciUtility::GetInstance()->input_format = argv[++i];
+        qaUtility::GetInstance()->format = argv[++i];
         continue;
       }
-      if (std::string(argv[i]) == "-input-format" && i == argc - 1)
+      if (std::string(argv[i]) == "-format" && i == argc - 1)
       {
         std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
         return 1;
       }
-      if (std::string(argv[i]) == "-output-format" && i != argc - 1)
+      if (std::string(argv[i]) == "-config" && i != argc - 1)
       {
-        mciUtility::GetInstance()->output_format = argv[++i];
+        configFileName = argv[++i];
         continue;
       }
-      if (std::string(argv[i]) == "-output-format" && i == argc - 1)
+      if (std::string(argv[i]) == "-config" && i == argc - 1)
       {
         std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
         return 1;
       }
-      if (std::string(argv[i]) == "-debug")
-      {
-        mciUtility::GetInstance()->debug = 1;
-        continue;
-      }
     }
   }
 
   TStopwatch timer;
   timer.Start();
-  
-  if (mciUtility::GetInstance()->debug)
+
+  Bool_t IsRead = true;
+
+  if (configFileName.Length() > 0)
   {
-    std::cout << "Input format = " << mciUtility::GetInstance()->input_format << std::endl;
-    std::cout << "Output format = " << mciUtility::GetInstance()->output_format << std::endl;
-    std::cout << "debug = " << mciUtility::GetInstance()->debug << std::endl;
-    std::cout << "Nevents = " << mciUtility::GetInstance()->Nevents << std::endl;
-    std::cout << std::endl;
+    IsRead = qaUtility::GetInstance()->ReadConfig(configFileName);
   }
 
-  mciReader_manager *readerManager = nullptr;
-  if (mciUtility::GetInstance()->input_format == "mcpico")
+  if (!IsRead)
   {
-    readerManager = new mciReader_mcpico();
+    std::cerr << "Error while reading config file. Exiting..." << std::endl;
+    return 2;
+  }
+
+  TFile *fo = new TFile(oFileName, "recreate");
+
+  const int ntracks_max = 15000;
+
+   Float_t         d_bimp;
+   Float_t         d_phi2;
+   Float_t         d_phi3;
+   Float_t         d_ecc2;
+   Float_t         d_ecc3;
+   Int_t           d_npart;
+   Int_t           d_nh;
+   Float_t         d_momx[ntracks_max];   //[nh]
+   Float_t         d_momy[ntracks_max];   //[nh]
+   Float_t         d_momz[ntracks_max];   //[nh]
+   Float_t         d_ene[ntracks_max];   //[nh]
+   Int_t           d_hid[ntracks_max];   //[nh]
+   Int_t           d_pdg[ntracks_max];   //[nh]
+   Short_t         d_charge[ntracks_max];   //[nh]
+
+  TTree *oTree = new TTree("mctree","mcpico format");
+  oTree->Branch("bimp",&d_bimp,"bimp/F"); // impact parametr
+  oTree->Branch("phi2",&d_phi2,"phi2/F"); // phiRP2
+  oTree->Branch("phi3",&d_phi3,"phi3/F"); // phiRP3
+  oTree->Branch("ecc2",&d_ecc2,"ecc2/F"); // eccentricity2
+  oTree->Branch("ecc3",&d_ecc3,"ecc3/F"); // eccentricity3
+  oTree->Branch("npart",&d_npart,"npart/I"); // number of participants
+  oTree->Branch("nh",&d_nh,"nh/I");  // number of particles
+  oTree->Branch("momx",&d_momx,"momx[nh]/F");
+  oTree->Branch("momy",&d_momy,"momy[nh]/F");
+  oTree->Branch("momz",&d_momz,"momz[nh]/F");
+  oTree->Branch("ene",&d_ene,"ene[nh]/F");//[energy]
+  oTree->Branch("hid",&d_hid,"hid[nh]/I");//[histrory id]
+  oTree->Branch("pdg",&d_pdg,"pdg[nh]/I");//[particle data group code]
+  oTree->Branch("charge",&d_charge,"charge[nh]/S");//[electric charge]
+
+  qaReader_manager *readerManager;
+  if (qaUtility::GetInstance()->format == "mcpico")
+  {
+    readerManager = new qaReader_mcpico();
   }
 #ifdef _MCINI_
-  if (mciUtility::GetInstance()->input_format == "mcini")
+  if (qaUtility::GetInstance()->format == "mcini")
   {
-    readerManager = new mciReader_mcini();
+    readerManager = new qaReader_mcini();
   }
 #endif
 #ifdef _PHQMD_
-  if (mciUtility::GetInstance()->input_format == "phqmd")
+  if (qaUtility::GetInstance()->format == "phqmd")
   {
-    readerManager = new mciReader_phqmd();
+    readerManager = new qaReader_phqmd();
   }
 #endif
 #ifdef _HSD_ROOT_
-  if (mciUtility::GetInstance()->input_format == "hsd")
+  if (qaUtility::GetInstance()->format == "hsd")
   {
-    readerManager = new mciReader_hsd_root();
+    readerManager = new qaReader_hsd_root();
   }
 #endif
-  if (mciUtility::GetInstance()->input_format == "particles")
+  if (qaUtility::GetInstance()->format == "particles")
   {
-    readerManager = new mciReader_smash_root();
+    readerManager = new qaReader_smash_root();
   }
 
   if (!readerManager)
@@ -157,54 +192,81 @@ int main(int argc, char **argv)
   readerManager->SetChain(iFileName.Data());
 
   Long64_t Nentries_chain = readerManager->GetEntries();
-  Long64_t Nentries = (mciUtility::GetInstance()->Nevents > Nentries_chain) ? Nentries_chain : mciUtility::GetInstance()->Nevents;
-  if (mciUtility::GetInstance()->Nevents == -1)
+  Long64_t Nentries = (qaUtility::GetInstance()->Nevents > Nentries_chain) ? Nentries_chain : qaUtility::GetInstance()->Nevents;
+  if (qaUtility::GetInstance()->Nevents == -1)
     Nentries = Nentries_chain;
   Int_t Nparticles;
 
-  mciWriter_manager *writerManager = nullptr;
-  if (mciUtility::GetInstance()->output_format == "mcpico")
-  {
-    writerManager = new mciWriter_mcpico();
-  }
-  if (!writerManager)
+  qaEvent *event = nullptr;
+  qaParticle *particle = nullptr;
+
+  Long64_t Absolute_counter = 0;
+  Long64_t Minbias_counter = 0;
+
+  while (Minbias_counter < Nentries)
   {
-    std::cerr << "This output format is not found!" << std::endl;
-    return 30;
-  }
+    //if (Minbias_counter % 1000 == 0)
+    std::cout << "Event [" << Minbias_counter << "/" << Nentries << "]" << std::endl;
 
-  writerManager->SetFile(oFileName.Data());
+    event = (qaEvent *)readerManager->ReadEvent(Absolute_counter);
+    Absolute_counter++;
 
-  mciEvent *event = nullptr;
-  mciParticle *particle = nullptr;
+    if (Absolute_counter > Nentries_chain)
+      break;
 
-  for (Long64_t ievent=0; ievent<Nentries; ievent++)
-  {
-    std::cout << "Event [" << ievent << "/" << Nentries << "]" << std::endl;
-    
-    event = (mciEvent *)readerManager->ReadEvent(ievent);
-    if (!event) continue;
+    if (!event)
+      continue;
+
+    Minbias_counter++;
 
-    writerManager->SetEvent(event);
+    d_bimp = event->GetB();
+    d_phi2 = event->GetPhiRP();
+    d_nh   = event->GetNparticles();
+    d_npart= 0;
+    d_ecc2 = 0.;
+    d_ecc3 = 0.;
 
     Nparticles = event->GetNparticles();
     for (int iparticle = 0; iparticle < Nparticles; iparticle++)
     {
-      particle = (mciParticle*) readerManager->ReadParticle(iparticle);
-      if (!particle) continue;
+      particle = readerManager->ReadParticle(iparticle);
 
-      writerManager->SetParticle(particle, iparticle);
-    }
+      if (!particle)
+      {
+        particle = new qaParticle();
+        particle->SetPx(qaUtility::GetInstance()->error_code);
+        particle->SetPy(qaUtility::GetInstance()->error_code);
+        particle->SetPz(qaUtility::GetInstance()->error_code);
+        particle->SetX(qaUtility::GetInstance()->error_code);
+        particle->SetY(qaUtility::GetInstance()->error_code);
+        particle->SetZ(qaUtility::GetInstance()->error_code);
+      }
 
-    writerManager->Fill();
+      d_momx[iparticle] = particle->GetPx();
+      d_momy[iparticle] = particle->GetPy();
+      d_momz[iparticle] = particle->GetPz();
+      d_ene[iparticle]  = particle->GetEnergy();  
+      d_pdg[iparticle]  = particle->GetPdg();  
+      d_hid[iparticle]  = 0.;
+      d_charge[iparticle] = particle->GetCharge();  
+
+      delete particle;
+    }
+    oTree->Fill();
+    delete event;
   }
 
-  writerManager->Write();
-  writerManager->Close();
+  std::cout << "Loop is closed, " << Minbias_counter << " events were counted." << std::endl;
+
+  std::cout << "Writing output." << std::endl;
+
+  fo->cd();
+  oTree->Write();
+  oTree->Print(); 
+  fo->Close();
 
   timer.Stop();
   timer.Print();
 
   return 0;
-
 }

format/qaEvent.cxx

@@ -0,0 +1,11 @@
+#include <qaEvent.h>
+
+ClassImp(qaEvent);
+
+qaEvent::qaEvent(/* args */)
+{
+}
+
+qaEvent::~qaEvent()
+{
+}

format/qaEvent.h

@@ -0,0 +1,33 @@
+#ifndef QATOOLS_FORMAT_EVENT_H
+#define QATOOLS_FORMAT_EVENT_H
+
+#include <iostream>
+#include <vector>
+
+#include <Rtypes.h>
+
+class qaEvent
+{
+private:
+  Float_t fB;
+  Float_t fPhiRP;
+  Int_t fNparticles;
+
+public:
+  qaEvent(/* args */);
+  virtual ~qaEvent();
+
+  // Setters
+  virtual void SetB(Float_t _a) { fB = _a; }
+  virtual void SetPhiRP(Float_t _a) { fPhiRP = _a; }
+  virtual void SetNparticles(Float_t _a) { fNparticles = _a; }
+
+  // Getters
+  virtual Float_t GetB() const { return fB; }
+  virtual Float_t GetPhiRP() const { return fPhiRP; }
+  virtual Int_t GetNparticles() const { return fNparticles; }
+
+  ClassDef(qaEvent, 0);
+};
+
+#endif

format/qaParticle.cxx

@@ -0,0 +1,11 @@
+#include <qaParticle.h>
+
+ClassImp(qaParticle);
+
+qaParticle::qaParticle(/* args */)
+{
+}
+
+qaParticle::~qaParticle()
+{
+}

format/qaParticle.h

@@ -0,0 +1,58 @@
+#ifndef QATOOLS_FORMAT_PARTICLE_H
+#define QATOOLS_FORMAT_PARTICLE_H
+
+#include <iostream>
+
+#include <Rtypes.h>
+#include <TVector3.h>
+
+class qaParticle
+{
+private:
+  TVector3 fMom;
+  TVector3 fR;
+  Float_t fEn;
+  Float_t fT;
+  Int_t fPdg;
+  Int_t fCharge;
+
+public:
+  qaParticle(/* args */);
+  virtual ~qaParticle();
+
+  // Setters
+  virtual void SetPxPyPz(Float_t _px, Float_t _py, Float_t _pz) { fMom.SetXYZ(_px, _py, _pz); }
+  virtual void SetXYZ(Float_t _x, Float_t _y, Float_t _z) { fR.SetXYZ(_x, _y, _z); }
+  virtual void SetPx(Float_t _a) { fMom.SetX(_a); }
+  virtual void SetPy(Float_t _a) { fMom.SetY(_a); }
+  virtual void SetPz(Float_t _a) { fMom.SetZ(_a); }
+  virtual void SetX(Float_t _a) { fR.SetX(_a); }
+  virtual void SetY(Float_t _a) { fR.SetY(_a); }
+  virtual void SetZ(Float_t _a) { fR.SetZ(_a); }
+  virtual void SetPdg(Int_t _a) { fPdg = _a; }
+  virtual void SetEnergy(Float_t _a) { fEn = _a; }
+  virtual void SetTime(Float_t _a) { fT = _a; }
+  virtual void SetCharge(Int_t _a) { fCharge = _a; }
+
+  // Getters
+  virtual TVector3 GetMom() { return fMom; }
+  virtual TVector3 GetR() { return fR; }
+  virtual Float_t GetPx() { return fMom.X(); }
+  virtual Float_t GetPy() { return fMom.Y(); }
+  virtual Float_t GetPz() { return fMom.Z(); }
+  virtual Float_t GetPt() { return fMom.Perp(); }
+  virtual Float_t GetP() { return fMom.Mag(); }
+  virtual Float_t GetEta() { return fMom.Eta(); }
+  virtual Float_t GetPhi() { return fMom.Phi(); }
+  virtual Float_t GetT() { return fT; }
+  virtual Float_t GetX() { return fR.X(); }
+  virtual Float_t GetY() { return fR.Y(); }
+  virtual Float_t GetZ() { return fR.Z(); }
+  virtual Int_t GetPdg() const { return fPdg; }
+  virtual Float_t GetEnergy() const { return fEn; }
+  virtual Int_t GetCharge() const { return fCharge; }
+
+  ClassDef(qaParticle, 0);
+};
+
+#endif

format/qaParticleLight.cxx

@@ -0,0 +1,31 @@
+#include <qaParticleLight.h>
+
+ClassImp(qaParticleLight);
+
+qaParticleLight::qaParticleLight(/* args */) : fPx(0), fPy(0), fPz(0), fE(0), fPdg(0)
+{
+}
+
+qaParticleLight::~qaParticleLight()
+{
+}
+
+void qaParticleLight::SetParticle(Float_t _px, Float_t _py, Float_t _pz, Float_t _e, Int_t _pdg, Int_t _charge)
+{
+    fPx = _px;
+    fPy = _py;
+    fPz = _pz;
+    fE = _e;
+    fPdg = _pdg;
+    fCharge = _charge;
+}
+
+void qaParticleLight::SetParticle(qaParticle *const &_particle)
+{
+    fPx = _particle->GetPx();
+    fPy = _particle->GetPy();
+    fPz = _particle->GetPz();
+    fE = _particle->GetEnergy();
+    fPdg = _particle->GetPdg();
+    fCharge = _particle->GetCharge();
+}

format/qaParticleLight.h

@@ -0,0 +1,43 @@
+#ifndef QATOOLS_FORMAT_PARTICLE_LIGHT_H
+#define QATOOLS_FORMAT_PARTICLE_LIGHT_H
+
+#include <iostream>
+#include <Rtypes.h>
+#include <TMath.h>
+
+#include <qaParticle.h>
+
+class qaParticleLight
+{
+private:
+    Float_t fPx;
+    Float_t fPy;
+    Float_t fPz;
+    Float_t fE;
+    Int_t fPdg;
+    Int_t fCharge;
+
+public:
+    qaParticleLight(/* args */);
+    virtual ~qaParticleLight();
+
+    virtual void SetParticle(Float_t _px, Float_t _py, Float_t _pz, Float_t _e, Int_t _pdg, Int_t _charge);
+    virtual void SetParticle(qaParticle *const &_particle);
+
+    virtual Float_t GetPx() const { return fPx; }
+    virtual Float_t GetPy() const { return fPy; }
+    virtual Float_t GetPz() const { return fPz; }
+    virtual Float_t GetEnergy() const { return fE; }
+    virtual Int_t GetPdg() const { return fPdg; }
+    virtual Int_t GetCharge() const { return fCharge; }
+
+    virtual Float_t GetPt() const { return TMath::Sqrt(TMath::Power(fPx, 2) + TMath::Power(fPy, 2)); }
+    virtual Float_t GetP() const { return TMath::Sqrt(TMath::Power(fPx, 2) + TMath::Power(fPy, 2) + TMath::Power(fPz, 2)); }
+    virtual Float_t GetEta() const { return 0.5 * TMath::Log((TMath::Sqrt(TMath::Power(fPx, 2) + TMath::Power(fPy, 2) + TMath::Power(fPz, 2)) + fPz) / (TMath::Sqrt(TMath::Power(fPx, 2) + TMath::Power(fPy, 2) + TMath::Power(fPz, 2)) - fPz)); }
+    virtual Float_t GetRapidity() const { return 0.5 * TMath::Log((fE + fPz) / (fE - fPz)); }
+    virtual Float_t GetPhi() const { return TMath::ATan2(fPy, fPx); }
+
+    ClassDef(qaParticleLight, 0);
+};
+
+#endif

qaTools.LinkDef.h

@@ -0,0 +1,28 @@
+#ifdef __CINT__
+
+#pragma link off all global;
+#pragma link off all classes;
+#pragma link off all functions;
+//#pragma link off all extern;
+
+#pragma link C++ nestedclass;
+#pragma link C++ nestedtypedef;
+
+#pragma link C++ class qaParticle+;
+#pragma link C++ class qaParticleLight+;
+#pragma link C++ class qaEvent+;
+#pragma link C++ class qaReader_smash_root+;
+#pragma link C++ class qaReader_mcpico+;
+#ifdef _MCINI_
+#pragma link C++ class qaReader_mcini+;
+#endif
+#ifdef _PHQMD_
+#pragma link C++ class qaReader_phqmd+;
+#endif
+#ifdef _HSD_ROOT_
+#pragma link C++ class qaReader_hsd_root+;
+#endif
+#pragma link C++ class qaReader_manager+;
+#pragma link C++ class qaUtility+;
+
+#endif

readers/qaReader_hsd_root.cxx

@@ -0,0 +1,136 @@
+#include <qaReader_hsd_root.h>
+
+ClassImp(qaReader_hsd_root);
+
+qaReader_hsd_root::qaReader_hsd_root(/* args */) : is_init(false), fCurrentEvent(-1), fEvent(nullptr), fParticle(nullptr)
+{
+}
+
+qaReader_hsd_root::~qaReader_hsd_root()
+{
+}
+
+Bool_t qaReader_hsd_root::ChainCheck()
+{
+  if (!is_init)
+  {
+    return false;
+  }
+
+  if (fCurrentEvent == -1)
+  {
+    return false;
+  }
+
+  if (!fChainPHSD->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+  if (!fChainFrag->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+
+  return true;
+}
+
+void qaReader_hsd_root::SetChain(const TString &inputFileName)
+{
+  fChainPHSD = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainPHSDName);
+
+  fChainPHSD->SetBranchAddress("event", &fEvent);
+
+  fChainFrag = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainFragName);
+
+  fChainFrag->SetBranchAddress("eventFrag", &fEventFrag);
+
+  is_init = kTRUE;
+}
+
+qaEvent *qaReader_hsd_root::ReadEvent(Long64_t iev)
+{
+  fCurrentEvent = iev;
+
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  qaEvent *event = new qaEvent();
+
+  event->SetB(fEvent->GetB());
+  event->SetPhiRP(0.); // in PHQMD axis orientation is rotated by 180 deg. In HSD/PHSD it's ok.
+  //event->SetPhiRP(TMath::Pi()); // in PHQMD axis orientation is rotated by 180 deg. In HSD/PHSD it's ok.
+  event->SetNparticles(fEvent->GetNpart() + fEventFrag->GetNfragMST()); // number of particles from PHSD and number of MST fragments
+
+  return event;
+}
+
+qaParticle *qaReader_hsd_root::ReadParticle(Int_t ipart)
+{
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  if (ipart >= fEvent->GetNpart() + fEventFrag->GetNfragMST())
+  {
+    return nullptr;
+  }
+
+  if (ipart < 0)
+  {
+    return nullptr;
+  }
+
+  bool is_particle = false;
+  bool is_fragment = false;
+
+  if (ipart < fEvent->GetNpart())
+  {
+    is_particle = true;
+  }
+  if (ipart >= fEvent->GetNpart())
+  {
+    is_fragment = true;
+  }
+
+  if (!is_particle && !is_fragment)
+  {
+    return nullptr;
+  }
+
+  if (is_particle && is_fragment)
+  {
+    return nullptr;
+  }
+
+  qaParticle *particle = new qaParticle();
+
+  if (is_particle)
+  {
+    fParticle = fEvent->GetParticle(ipart);
+    if (fParticle->IsInMST()) // Check whether the particle from PHSD was used to create MST fragment
+    {
+      return nullptr;
+    }
+    particle->SetEnergy(fParticle->E());
+    particle->SetPdg(fParticle->GetPdg());
+    particle->SetPxPyPz(fParticle->Px(), fParticle->Py(), fParticle->Pz());
+    particle->SetTime(0.);
+    particle->SetXYZ(0., 0., 0.);
+    particle->SetCharge(fParticle->GetCharge());
+  }
+  if (is_fragment)
+  {
+    fFragmentMST = fEventFrag->GetFragmentMST(ipart - fEvent->GetNpart());
+    particle->SetEnergy(fFragmentMST->E());
+    particle->SetPdg(fFragmentMST->GetPdg());
+    particle->SetPxPyPz(fFragmentMST->Px(), fFragmentMST->Py(), fFragmentMST->Pz());
+    particle->SetTime(0.);
+    particle->SetXYZ(0., 0., 0.);
+    particle->SetCharge(fFragmentMST->GetZ());
+  }
+
+  return particle;
+}

readers/qaReader_hsd_root.h

@@ -0,0 +1,48 @@
+#ifndef QATOOLS_READERS_HSD_ROOT_H
+#define QATOOLS_READERS_HSD_ROOT_H
+
+#include <Rtypes.h>
+#include <TChain.h>
+#include <TLorentzVector.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+#include <qaReader_manager.h>
+
+#include <phqmd2root/src/libPHQMDEvent.h>
+
+class qaReader_hsd_root : virtual public qaReader_manager
+{
+private:
+  TChain *fChainPHSD, *fChainFrag;
+
+  const char *fChainPHSDName = "PHQMDtree";
+  const char *fChainFragName = "FRAGtree";
+
+  bool is_init;
+  Long64_t fCurrentEvent;
+
+  Event *fEvent;
+  EventClust *fEventFrag;
+  Particle *fParticle;
+  Baryon *fBaryonMST;
+  Fragment *fFragmentMST;
+  TLorentzVector fMomentum;
+
+  Bool_t ChainCheck();
+
+public:
+  qaReader_hsd_root(/* args */);
+  virtual ~qaReader_hsd_root();
+
+  virtual void SetChain(const TString &inputFileName);
+  virtual Long64_t GetEntries() { return fChainFrag->GetEntries(); }
+  virtual qaEvent *ReadEvent(Long64_t iev);
+  virtual qaParticle* ReadParticle(Int_t ipart);
+
+  ClassDef(qaReader_hsd_root, 0);
+};
+
+#endif

readers/qaReader_manager.cxx

@@ -0,0 +1,7 @@
+#include <qaReader_manager.h>
+
+ClassImp(qaReader_manager);
+
+qaReader_manager::~qaReader_manager()
+{
+}

readers/qaReader_manager.h

@@ -0,0 +1,26 @@
+#ifndef QATOOLS_READERS_MANAGER_H
+#define QATOOLS_READERS_MANAGER_H
+
+#include <Rtypes.h>
+#include <TString.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+class qaReader_manager
+{
+private:
+  /* data */
+public:
+  virtual ~qaReader_manager();
+
+  virtual void SetChain(const TString &inputFileName) = 0;
+  virtual Long64_t GetEntries() = 0;
+  virtual qaEvent *ReadEvent(Long64_t iev) = 0;
+  virtual qaParticle* ReadParticle(Int_t ipart) = 0;
+
+  ClassDef(qaReader_manager,0);
+};
+
+#endif

readers/qaReader_mcini.cxx

@@ -0,0 +1,85 @@
+#include <qaReader_mcini.h>
+
+ClassImp(qaReader_mcini);
+
+qaReader_mcini::qaReader_mcini(/* args */) : is_init(false), fCurrentEvent(-1), fEvent(nullptr), fParticle(nullptr)
+{
+}
+
+qaReader_mcini::~qaReader_mcini()
+{
+}
+
+Bool_t qaReader_mcini::ChainCheck()
+{
+  if (!is_init)
+  {
+    return false;
+  }
+
+  if (fCurrentEvent == -1)
+  {
+    return false;
+  }
+
+  if (!fChain->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+
+  return true;
+}
+
+void qaReader_mcini::SetChain(const TString &inputFileName)
+{
+  fChain = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainName);
+
+  fChain->SetBranchAddress("event", &fEvent);
+
+  is_init = kTRUE;
+}
+
+qaEvent *qaReader_mcini::ReadEvent(Long64_t iev)
+{
+  fCurrentEvent = iev;
+
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  qaEvent *event = new qaEvent();
+
+  event->SetB(fEvent->GetB());
+  event->SetPhiRP(fEvent->GetPhi());
+  event->SetNparticles(fEvent->GetNpa());
+
+  return event;
+}
+
+qaParticle *qaReader_mcini::ReadParticle(Int_t ipart)
+{
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  if (ipart >= fEvent->GetNpa())
+  {
+    return nullptr;
+  }
+
+  fParticle = fEvent->GetParticle(ipart);
+  fMomentum = fParticle->GetMomentum();
+
+  qaParticle *particle = new qaParticle();
+
+  particle->SetEnergy(fMomentum.E());
+  particle->SetPdg(fParticle->GetPdg());
+  particle->SetPxPyPz(fMomentum.Px(), fMomentum.Py(), fMomentum.Pz());
+  particle->SetTime(0.);
+  particle->SetXYZ(0., 0., 0.);
+  particle->SetCharge(qaUtility::GetInstance()->GetCharge(fParticle->GetPdg()));
+
+  return particle;
+}

readers/qaReader_mcini.h

@@ -0,0 +1,47 @@
+#ifndef QATOOLS_READERS_MCINI_H
+#define QATOOLS_READERS_MCINI_H
+
+#include <Rtypes.h>
+#include <TChain.h>
+#include <TLorentzVector.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+#include <qaReader_manager.h>
+
+#include <URun.h>
+#include <UEvent.h>
+#include <UParticle.h>
+#include <EventInitialState.h>
+
+class qaReader_mcini : virtual public qaReader_manager
+{
+private:
+  TChain *fChain;
+
+  const char *fChainName = "events";
+
+  bool is_init;
+  Long64_t fCurrentEvent;
+
+  UEvent *fEvent;
+  UParticle *fParticle;
+  TLorentzVector fMomentum;
+
+  Bool_t ChainCheck();
+
+public:
+  qaReader_mcini(/* args */);
+  virtual ~qaReader_mcini();
+
+  virtual void SetChain(const TString &inputFileName);
+  virtual Long64_t GetEntries() { return fChain->GetEntries(); }
+  virtual qaEvent *ReadEvent(Long64_t iev);
+  virtual qaParticle* ReadParticle(Int_t ipart);
+
+  ClassDef(qaReader_mcini, 0);
+};
+
+#endif

readers/qaReader_mcpico.cxx

@@ -0,0 +1,95 @@
+#include <qaReader_mcpico.h>
+
+ClassImp(qaReader_mcpico);
+
+qaReader_mcpico::qaReader_mcpico(/* args */) : is_init(false), fCurrentEvent(-1)
+{
+}
+
+qaReader_mcpico::~qaReader_mcpico()
+{
+}
+
+Bool_t qaReader_mcpico::ChainCheck()
+{
+  if (!is_init)
+  {
+    return false;
+  }
+
+  if (fCurrentEvent == -1)
+  {
+    return false;
+  }
+
+  if (!fChain->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+
+  return true;
+}
+
+void qaReader_mcpico::SetChain(const TString &inputFileName)
+{
+  fChain = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainName);
+
+  fChain->SetBranchAddress("bimp", &bimp);
+  fChain->SetBranchAddress("phi2", &phi2);
+  fChain->SetBranchAddress("phi3", &phi3);
+  fChain->SetBranchAddress("ecc2", &ecc2);
+  fChain->SetBranchAddress("ecc3", &ecc3);
+  fChain->SetBranchAddress("npart", &npart);
+  fChain->SetBranchAddress("nh", &nh);
+  fChain->SetBranchAddress("momx", momx);
+  fChain->SetBranchAddress("momy", momy);
+  fChain->SetBranchAddress("momz", momz);
+  fChain->SetBranchAddress("ene", ene);
+  fChain->SetBranchAddress("hid", hid);
+  fChain->SetBranchAddress("pdg", pdg);
+  fChain->SetBranchAddress("charge", charge);
+
+  is_init = kTRUE;
+}
+
+qaEvent *qaReader_mcpico::ReadEvent(Long64_t iev)
+{
+  fCurrentEvent = iev;
+
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  qaEvent *event = new qaEvent();
+
+  event->SetB(bimp);
+  event->SetPhiRP(phi2);
+  event->SetNparticles(nh);
+
+  return event;
+}
+
+qaParticle *qaReader_mcpico::ReadParticle(Int_t ipart)
+{
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  if (ipart >= nh)
+  {
+    return nullptr;
+  }
+
+  qaParticle *particle = new qaParticle();
+
+  particle->SetEnergy(ene[ipart]);
+  particle->SetPdg(pdg[ipart]);
+  particle->SetPxPyPz(momx[ipart], momy[ipart], momz[ipart]);
+  particle->SetTime(0.);
+  particle->SetXYZ(0., 0., 0.);
+  particle->SetCharge(qaUtility::GetInstance()->GetCharge(pdg[ipart]));
+
+  return particle;
+}

readers/qaReader_mcpico.h

@@ -0,0 +1,54 @@
+#ifndef QATOOLS_READERS_MCPICO_H
+#define QATOOLS_READERS_MCPICO_H
+
+#include <Rtypes.h>
+#include <TChain.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+#include <qaReader_manager.h>
+
+#define MAX_TRACKS 15000
+
+class qaReader_mcpico : virtual public qaReader_manager
+{
+private:
+  TChain *fChain;
+
+  const char *fChainName = "mctree";
+
+  bool is_init;
+  Long64_t fCurrentEvent;
+
+   Float_t         bimp;
+   Float_t         phi2;
+   Float_t         phi3;
+   Float_t         ecc2;
+   Float_t         ecc3;
+   Int_t           npart;
+   Int_t           nh;
+   Float_t         momx[MAX_TRACKS];   //[nh]
+   Float_t         momy[MAX_TRACKS];   //[nh]
+   Float_t         momz[MAX_TRACKS];   //[nh]
+   Float_t         ene[MAX_TRACKS];   //[nh]
+   Int_t           hid[MAX_TRACKS];   //[nh]
+   Int_t           pdg[MAX_TRACKS];   //[nh]
+   Short_t         charge[MAX_TRACKS];   //[nh]
+
+  Bool_t ChainCheck();
+
+public:
+  qaReader_mcpico(/* args */);
+  virtual ~qaReader_mcpico();
+
+  virtual void SetChain(const TString &inputFileName);
+  virtual Long64_t GetEntries() { return fChain->GetEntries(); }
+  virtual qaEvent *ReadEvent(Long64_t iev);
+  virtual qaParticle* ReadParticle(Int_t ipart);
+
+  ClassDef(qaReader_mcpico, 0);
+};
+
+#endif

readers/qaReader_phqmd.cxx

@@ -0,0 +1,136 @@
+#include <qaReader_phqmd.h>
+
+ClassImp(qaReader_phqmd);
+
+qaReader_phqmd::qaReader_phqmd(/* args */) : is_init(false), fCurrentEvent(-1), fEvent(nullptr), fParticle(nullptr)
+{
+}
+
+qaReader_phqmd::~qaReader_phqmd()
+{
+}
+
+Bool_t qaReader_phqmd::ChainCheck()
+{
+  if (!is_init)
+  {
+    return false;
+  }
+
+  if (fCurrentEvent == -1)
+  {
+    return false;
+  }
+
+  if (!fChainPHSD->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+  if (!fChainFrag->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+
+  return true;
+}
+
+void qaReader_phqmd::SetChain(const TString &inputFileName)
+{
+  fChainPHSD = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainPHSDName);
+
+  fChainPHSD->SetBranchAddress("event", &fEvent);
+
+  fChainFrag = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainFragName);
+
+  fChainFrag->SetBranchAddress("eventFrag", &fEventFrag);
+
+  is_init = kTRUE;
+}
+
+qaEvent *qaReader_phqmd::ReadEvent(Long64_t iev)
+{
+  fCurrentEvent = iev;
+
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  qaEvent *event = new qaEvent();
+
+  event->SetB(fEvent->GetB());
+  // event->SetPhiRP(0.); // in PHQMD axis orientation is rotated by 180 deg. In HSD/PHSD it's ok.
+  event->SetPhiRP(TMath::Pi()); // in PHQMD axis orientation is rotated by 180 deg. In HSD/PHSD it's ok.
+  event->SetNparticles(fEvent->GetNpart() + fEventFrag->GetNfragMST()); // number of particles from PHSD and number of MST fragments
+
+  return event;
+}
+
+qaParticle *qaReader_phqmd::ReadParticle(Int_t ipart)
+{
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  if (ipart >= fEvent->GetNpart() + fEventFrag->GetNfragMST())
+  {
+    return nullptr;
+  }
+
+  if (ipart < 0)
+  {
+    return nullptr;
+  }
+
+  bool is_particle = false;
+  bool is_fragment = false;
+
+  if (ipart < fEvent->GetNpart())
+  {
+    is_particle = true;
+  }
+  if (ipart >= fEvent->GetNpart())
+  {
+    is_fragment = true;
+  }
+
+  if (!is_particle && !is_fragment)
+  {
+    return nullptr;
+  }
+
+  if (is_particle && is_fragment)
+  {
+    return nullptr;
+  }
+
+  qaParticle *particle = new qaParticle();
+
+  if (is_particle)
+  {
+    fParticle = fEvent->GetParticle(ipart);
+    if (fParticle->IsInMST()) // Check whether the particle from PHSD was used to create MST fragment
+    {
+      return nullptr;
+    }
+    particle->SetEnergy(fParticle->E());
+    particle->SetPdg(fParticle->GetPdg());
+    particle->SetPxPyPz(fParticle->Px(), fParticle->Py(), fParticle->Pz());
+    particle->SetTime(0.);
+    particle->SetXYZ(0., 0., 0.);
+    particle->SetCharge(fParticle->GetCharge());
+  }
+  if (is_fragment)
+  {
+    fFragmentMST = fEventFrag->GetFragmentMST(ipart - fEvent->GetNpart());
+    particle->SetEnergy(fFragmentMST->E());
+    particle->SetPdg(fFragmentMST->GetPdg());
+    particle->SetPxPyPz(fFragmentMST->Px(), fFragmentMST->Py(), fFragmentMST->Pz());
+    particle->SetTime(0.);
+    particle->SetXYZ(0., 0., 0.);
+    particle->SetCharge(fFragmentMST->GetZ());
+  }
+
+  return particle;
+}

readers/qaReader_phqmd.h

@@ -0,0 +1,48 @@
+#ifndef QATOOLS_READERS_PHQMD_H
+#define QATOOLS_READERS_PHQMD_H
+
+#include <Rtypes.h>
+#include <TChain.h>
+#include <TLorentzVector.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+#include <qaReader_manager.h>
+
+#include <phqmd2root/src/libPHQMDEvent.h>
+
+class qaReader_phqmd : virtual public qaReader_manager
+{
+private:
+  TChain *fChainPHSD, *fChainFrag;
+
+  const char *fChainPHSDName = "PHQMDtree";
+  const char *fChainFragName = "FRAGtree";
+
+  bool is_init;
+  Long64_t fCurrentEvent;
+
+  Event *fEvent;
+  EventClust *fEventFrag;
+  Particle *fParticle;
+  Baryon *fBaryonMST;
+  Fragment *fFragmentMST;
+  TLorentzVector fMomentum;
+
+  Bool_t ChainCheck();
+
+public:
+  qaReader_phqmd(/* args */);
+  virtual ~qaReader_phqmd();
+
+  virtual void SetChain(const TString &inputFileName);
+  virtual Long64_t GetEntries() { return fChainFrag->GetEntries(); }
+  virtual qaEvent *ReadEvent(Long64_t iev);
+  virtual qaParticle* ReadParticle(Int_t ipart);
+
+  ClassDef(qaReader_phqmd, 0);
+};
+
+#endif

readers/qaReader_smash_root.cxx

@@ -0,0 +1,92 @@
+#include <qaReader_smash_root.h>
+
+ClassImp(qaReader_smash_root);
+
+qaReader_smash_root::qaReader_smash_root(/* args */) : is_init(false), fCurrentEvent(-1)
+{
+}
+
+qaReader_smash_root::~qaReader_smash_root()
+{
+}
+
+Bool_t qaReader_smash_root::ChainCheck()
+{
+  if (!is_init)
+  {
+    return false;
+  }
+
+  if (fCurrentEvent == -1)
+  {
+    return false;
+  }
+
+  if (!fChain->GetEntry(fCurrentEvent))
+  {
+    return false;
+  }
+
+  return true;
+}
+
+void qaReader_smash_root::SetChain(const TString &inputFileName)
+{
+  fChain = (TChain *)qaUtility::GetInstance()->initChain(inputFileName, fChainName);
+
+  fChain->SetBranchAddress("npart", &npart);
+  fChain->SetBranchAddress("impact_b", &impact_b);
+  fChain->SetBranchAddress("pdgcode", pdgcode);
+  fChain->SetBranchAddress("p0", p0);
+  fChain->SetBranchAddress("px", px);
+  fChain->SetBranchAddress("py", py);
+  fChain->SetBranchAddress("pz", pz);
+  fChain->SetBranchAddress("t", t);
+  fChain->SetBranchAddress("x", x);
+  fChain->SetBranchAddress("y", y);
+  fChain->SetBranchAddress("z", z);
+
+  is_init = kTRUE;
+}
+
+qaEvent *qaReader_smash_root::ReadEvent(Long64_t iev)
+{
+  fCurrentEvent = iev;
+
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  qaEvent *event = new qaEvent();
+
+  event->SetB(impact_b);
+  event->SetPhiRP(0.);
+  event->SetNparticles(npart);
+
+  return event;
+}
+
+qaParticle *qaReader_smash_root::ReadParticle(Int_t ipart)
+{
+  if (!ChainCheck())
+  {
+    return nullptr;
+  }
+
+  if (ipart >= npart)
+  {
+    return nullptr;
+  }
+
+  qaParticle *particle = new qaParticle();
+
+  particle->SetEnergy(p0[ipart]);
+  particle->SetPdg(pdgcode[ipart]);
+  particle->SetPxPyPz(px[ipart], py[ipart], pz[ipart]);
+  particle->SetTime(t[ipart]);
+  particle->SetXYZ(x[ipart], y[ipart], z[ipart]);
+  particle->SetCharge(qaUtility::GetInstance()->GetCharge(pdgcode[ipart]));
+
+  return particle;
+}

readers/qaReader_smash_root.h

@@ -0,0 +1,53 @@
+#ifndef QATOOLS_READERS_SMASH_ROOT_H
+#define QATOOLS_READERS_SMASH_ROOT_H
+
+#include <Rtypes.h>
+#include <TChain.h>
+
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <Utility.h>
+
+#include <qaReader_manager.h>
+
+#define MAX_TRACKS 15000
+
+class qaReader_smash_root : virtual public qaReader_manager
+{
+private:
+  TChain *fChain;
+
+  const char *fChainName = "particles";
+
+  bool is_init;
+  Long64_t fCurrentEvent;
+
+  Int_t npart;
+  Double_t impact_b;
+  Int_t ev;
+  Int_t tcounter;
+  Int_t pdgcode[MAX_TRACKS]; //[npart]
+  Double_t p0[MAX_TRACKS];   //[npart]
+  Double_t px[MAX_TRACKS];   //[npart]
+  Double_t py[MAX_TRACKS];   //[npart]
+  Double_t pz[MAX_TRACKS];   //[npart]
+  Double_t t[MAX_TRACKS];    //[npart]
+  Double_t x[MAX_TRACKS];    //[npart]
+  Double_t y[MAX_TRACKS];    //[npart]
+  Double_t z[MAX_TRACKS];    //[npart]
+
+  Bool_t ChainCheck();
+
+public:
+  qaReader_smash_root(/* args */);
+  virtual ~qaReader_smash_root();
+
+  virtual void SetChain(const TString &inputFileName);
+  virtual Long64_t GetEntries() { return fChain->GetEntries(); }
+  virtual qaEvent *ReadEvent(Long64_t iev);
+  virtual qaParticle* ReadParticle(Int_t ipart);
+
+  ClassDef(qaReader_smash_root, 0);
+};
+
+#endif

utility/Utility.cxx

@@ -2,32 +2,52 @@
 #include <TDatabasePDG.h>
 #include <TParticlePDG.h>
 
-ClassImp(mciUtility);
+ClassImp(qaUtility);
 
-mciUtility *mciUtility::fUtility = nullptr;
+qaUtility *qaUtility::fUtility = nullptr;
 ;
 
-mciUtility::mciUtility() : Nevents(-1),
+qaUtility::qaUtility() : Nevents(-1),
                          debug(0),
-                         input_format("mcpico"),
-                         output_format("mcpico")
+                         format("mctree"),
+                         Is_minbias(1),
+                         Is_refmult(1),
+                         Is_v1(1),
+                         Is_v2(1),
+                         Is_v3(1),
+                         Is_v4(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_refmult_Particle_isCharged(1)
 {
 }
 
-mciUtility::~mciUtility()
+qaUtility::~qaUtility()
 {
 }
 
-mciUtility *mciUtility::GetInstance()
+qaUtility *qaUtility::GetInstance()
 {
   if (fUtility == nullptr)
   {
-    fUtility = new mciUtility();
+    fUtility = new qaUtility();
   }
   return fUtility;
 }
 
-Bool_t mciUtility::ReadConfig(const TString &configFileName)
+Bool_t qaUtility::ReadConfig(const TString &configFileName)
 {
   if (configFileName.Length() == 0)
   {
@@ -39,34 +59,49 @@ Bool_t mciUtility::ReadConfig(const TString &configFileName)
   Nevents = env.GetValue("Nevents", 0);
   debug = env.GetValue("debug", 0);
 
+  Is_minbias = env.GetValue("Is_minbias", 0);
+  Cut_minbias_Event_bmin = env.GetValue("Cut_minbias_Event_bmin", 0.);
+  Cut_minbias_Event_bmax = env.GetValue("Cut_minbias_Event_bmax", 0.);
+
+  Cut_minbias_Particle_ptmin = env.GetValue("Cut_minbias_Particle_ptmin", 0.);
+  Cut_minbias_Particle_ptmax = env.GetValue("Cut_minbias_Particle_ptmax", 0.);
+  Cut_minbias_Particle_etamin = env.GetValue("Cut_minbias_Particle_etamin", 0.);
+  Cut_minbias_Particle_etamax = env.GetValue("Cut_minbias_Particle_etamax", 0.);
+  Cut_minbias_Particle_ymin = env.GetValue("Cut_minbias_Particle_ymin", 0.);
+  Cut_minbias_Particle_ymax = env.GetValue("Cut_minbias_Particle_ymax", 0.);
+
+  Is_refmult = env.GetValue("Is_refmult", 0);
+  Cut_refmult_Event_bmin = env.GetValue("Cut_refmult_Event_bmin", 0.);
+  Cut_refmult_Event_bmax = env.GetValue("Cut_refmult_Event_bmax", 0.);
+
+  Cut_refmult_Particle_ptmin = env.GetValue("Cut_refmult_Particle_ptmin", 0.);
+  Cut_refmult_Particle_ptmax = env.GetValue("Cut_refmult_Particle_ptmax", 0.);
+  Cut_refmult_Particle_etamin = env.GetValue("Cut_refmult_Particle_etamin", 0.);
+  Cut_refmult_Particle_etamax = env.GetValue("Cut_refmult_Particle_etamax", 0.);
+  Cut_refmult_Particle_isCharged = env.GetValue("Cut_refmult_Particle_isCharged", 0);
+
+  Is_v1 = env.GetValue("Is_v1", 0);
+  Is_v2 = env.GetValue("Is_v2", 0);
+  Is_v3 = env.GetValue("Is_v3", 0);
+  Is_v4 = env.GetValue("Is_v4", 0);
+
   return true;
 }
 
-TChain *mciUtility::initChain(const TString &inputFileName, const char *chainName)
+TChain *qaUtility::initChain(const TString &inputFileName, const char *chainName)
 {
   TChain *chain = new TChain(chainName);
   std::ifstream file(inputFileName.Data());
   std::string line;
-  std::string basename = inputFileName.Data();
-  if (basename.find(".root") != std::string::npos)
-  {
-    chain->Add(inputFileName.Data());
-  }
-  if (file.is_open())
+  while (std::getline(file, line))
   {
-    while (std::getline(file, line))
-    {
-      if (line.find(".root") != std::string::npos)
-      {
-        chain->Add(line.c_str());
-      }
-    }
+    chain->Add(line.c_str());
   }
 
   return chain;
 }
 
-std::vector<Float_t> mciUtility::ParseVector(std::string _input)
+std::vector<Float_t> qaUtility::ParseVector(std::string _input)
 {
   std::vector<Float_t> vB;
 
@@ -79,10 +114,116 @@ std::vector<Float_t> mciUtility::ParseVector(std::string _input)
   return vB;
 }
 
-Double_t mciUtility::GetCharge(Int_t pdg)
+Bool_t qaUtility::Cut_Event_minbias(qaEvent *const &event)
+{
+  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_refmult(qaEvent *const &event)
+{
+  if (event->GetB() < Cut_refmult_Event_bmin)
+    return false;
+  if (event->GetB() > Cut_refmult_Event_bmax)
+    return false;
+
+  return true;
+}
+
+Bool_t qaUtility::Cut_Particle_minbias(qaParticle *const &particle)
+{
+  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 (Cut_refmult_Particle_isCharged)
+  {
+    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.;
 }
+
+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(4))
+    return 4;
+  if (pdg == vpdg.at(6))
+    return 6;
+  if (pdg == vpdg.at(7))
+    return 7;
+  if (pdg == vpdg.at(8))
+    return 8;
+  if (pdg == vpdg.at(9))
+    return 9;
+  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;
+}
+
+Int_t qaUtility::GetCentMultBin(Int_t mult, std::vector<Float_t> vcent)
+{
+  if (vcent.size() == 0)
+    return -1;
+
+  for (int i = 0; i < (int)vcent.size() - 1; i++)
+  {
+    if (mult >= vcent.at(i) && mult < vcent.at(i + 1))
+      return i;
+  }
+
+  return -1;
+}

utility/Utility.h

@@ -1,5 +1,5 @@
-#ifndef MODELCONVERTER_UTILITY_H
-#define MODELCONVERTER_UTILITY_H
+#ifndef QATOOLS_UTILITY_H
+#define QATOOLS_UTILITY_H
 
 #include <iostream>
 #include <fstream>
@@ -13,30 +13,58 @@
 #include <Rtypes.h>
 #include <TString.h>
 #include <TChain.h>
-#include <TFile.h>
 #include <TEnv.h>
 
-#define MAX_TRACKS 15000
+#include <qaEvent.h>
+#include <qaParticle.h>
+#include <qaParticleLight.h>
 
-class mciUtility
+class qaUtility
 {
 protected:
-  mciUtility();
-  virtual ~mciUtility();
-  static mciUtility *fUtility;
+  qaUtility();
+  virtual ~qaUtility();
+  static qaUtility *fUtility;
 
 public:
-  mciUtility(mciUtility &other) = delete;
-  void operator=(const mciUtility &) = delete;
+  qaUtility(qaUtility &other) = delete;
+  void operator=(const qaUtility &) = delete;
 
-  static mciUtility *GetInstance();
+  static qaUtility *GetInstance();
 
   const Double_t error_code = -999.;
+  const Int_t npid = 10;
+  const std::vector<Int_t> vpdg = {0, 211, 321, 2212, 2112, 0, -211, -321, -2212, -2112};
+  const std::vector<Double_t> mpdg = {error_code, 0.13957, 0.493677, 0.938272, 0.88277299, error_code, 0.13957, 0.493677, 0.938272, 0.88277299};
+  const Int_t maxCentBins = 20;
 
   Int_t Nevents;
   Int_t debug;
-  std::string input_format;
-  std::string output_format;
+  std::string format;
+
+  Int_t Is_minbias;
+  Int_t Is_refmult;
+  Int_t Is_v1;
+  Int_t Is_v2;
+  Int_t Is_v3;
+  Int_t Is_v4;
+
+  Double_t Cut_minbias_Event_bmin;
+  Double_t Cut_minbias_Event_bmax;
+  Double_t Cut_minbias_Particle_ptmin;
+  Double_t Cut_minbias_Particle_ptmax;
+  Double_t Cut_minbias_Particle_etamin;
+  Double_t Cut_minbias_Particle_etamax;
+  Double_t Cut_minbias_Particle_ymin;
+  Double_t Cut_minbias_Particle_ymax;
+
+  Double_t Cut_refmult_Event_bmin;
+  Double_t Cut_refmult_Event_bmax;
+  Double_t Cut_refmult_Particle_ptmin;
+  Double_t Cut_refmult_Particle_ptmax;
+  Double_t Cut_refmult_Particle_etamin;
+  Double_t Cut_refmult_Particle_etamax;
+  Int_t    Cut_refmult_Particle_isCharged;
 
   Bool_t ReadConfig(const TString &configFileName);
   TChain *initChain(const TString &inputFileName, const char *chainName);
@@ -44,9 +72,18 @@ public:
   std::vector<Float_t> ParseVector(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_Particle_minbias(qaParticle *const &particle);
+  Bool_t Cut_Particle_refmult(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);
+  Int_t GetCentMultBin(Int_t mult, std::vector<Float_t> vcent);
 
-  ClassDef(mciUtility, 0);
-}; // class mciUtility
+  ClassDef(qaUtility, 0);
+}; // class qaUtility
 
 #endif