ModelConverter/bin/main.cpp

#include <iostream>
#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 <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 <qaReader_mcini.h>
#endif
#ifdef _PHQMD_
#include <qaReader_phqmd.h>
#endif
#ifdef _HSD_ROOT_
#include <qaReader_hsd_root.h>
#endif

int main(int argc, char **argv)
{
  TString iFileName, oFileName, configFileName = "";

  if (argc < 7)
  {
    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_
    std::cerr << "\tmcini    - custom ROOT format to store both initial state and final state (UniGen data format) model data. Fully compatible with UniGen format." << std::endl;
  #endif
  #ifdef _PHQMD_
    std::cerr << "\tphqmd    - custom ROOT format to store PHQMD (with MST) model data." << std::endl;
  #endif
  #ifdef _HSD_ROOT_
    std::cerr << "\thsd      - custom ROOT format to store HSD model data." << std::endl;
  #endif
    return 1;
  }
  for (int i = 1; i < argc; i++)
  {
    if (std::string(argv[i]) != "-i" &&
        std::string(argv[i]) != "-o" &&
        std::string(argv[i]) != "-format" &&
        std::string(argv[i]) != "-config")
    {
      std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
      return 2;
    }
    else
    {
      if (std::string(argv[i]) == "-i" && i != argc - 1)
      {
        iFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-i" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "-o" && i != argc - 1)
      {
        oFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-o" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "-format" && i != argc - 1)
      {
        qaUtility::GetInstance()->format = argv[++i];
        continue;
      }
      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]) == "-config" && i != argc - 1)
      {
        configFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-config" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
    }
  }

  TStopwatch timer;
  timer.Start();

  Bool_t IsRead = true;

  if (configFileName.Length() > 0)
  {
    IsRead = qaUtility::GetInstance()->ReadConfig(configFileName);
  }

  if (!IsRead)
  {
    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 (qaUtility::GetInstance()->format == "mcini")
  {
    readerManager = new qaReader_mcini();
  }
#endif
#ifdef _PHQMD_
  if (qaUtility::GetInstance()->format == "phqmd")
  {
    readerManager = new qaReader_phqmd();
  }
#endif
#ifdef _HSD_ROOT_
  if (qaUtility::GetInstance()->format == "hsd")
  {
    readerManager = new qaReader_hsd_root();
  }
#endif
  if (qaUtility::GetInstance()->format == "particles")
  {
    readerManager = new qaReader_smash_root();
  }

  if (!readerManager)
  {
    std::cerr << "This input format is not found!" << std::endl;
    return 20;
  }

  readerManager->SetChain(iFileName.Data());

  Long64_t Nentries_chain = readerManager->GetEntries();
  Long64_t Nentries = (qaUtility::GetInstance()->Nevents > Nentries_chain) ? Nentries_chain : qaUtility::GetInstance()->Nevents;
  if (qaUtility::GetInstance()->Nevents == -1)
    Nentries = Nentries_chain;
  Int_t Nparticles;

  qaEvent *event = nullptr;
  qaParticle *particle = nullptr;

  Long64_t Absolute_counter = 0;
  Long64_t Minbias_counter = 0;

  while (Minbias_counter < Nentries)
  {
    //if (Minbias_counter % 1000 == 0)
    std::cout << "Event [" << Minbias_counter << "/" << Nentries << "]" << std::endl;

    event = (qaEvent *)readerManager->ReadEvent(Absolute_counter);
    Absolute_counter++;

    if (Absolute_counter > Nentries_chain)
      break;

    if (!event)
      continue;

    Minbias_counter++;

    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 = readerManager->ReadParticle(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);
      }

      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;
  }

  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;
}