KVSimReader_IQMD.cpp

#include "KVSimReader_IQMD.h"
#include "KV2Body.h"
 
 
 
KVSimReader_IQMD::KVSimReader_IQMD()
{
   SetTreeName("IQMD");
}
 
 
 
 
KVSimReader_IQMD::KVSimReader_IQMD(const TString& filename)
   : KVSimReader_IQMD{}
{
   ConvertEventsInFile(filename);
}
 
 
 
 
Bool_t KVSimReader_IQMD::ReadFileHeader()
{
   auto status = ReadLineAndCheck(6," "); //  30   40   83  126   34.600     zp,np,zt,nt,ein(MeV/n)
   if(status == ReadStatus::OK)
   {
      auto zp = GetIntReadPar(0);
      auto np = GetIntReadPar(1);
      auto zt = GetIntReadPar(2);
      auto nt = GetIntReadPar(3);
      auto ebeam = GetDoubleReadPar(4);
      KV2Body reac({zp,zp+np,(zp+np)*ebeam},{zt,zt+nt});
      KVError::Info(this,"ReadHeader", "Reading events for reaction: %s", reac.GetTitle());
 
      SetROOTFileName(Form("IQMD_%s_%s_%s_%.1fAMeV",
                           sim_type.Data(),
                           reac.GetNucleus(KV2Body::projectile)->GetSymbol(),
                           reac.GetNucleus(KV2Body::target)->GetSymbol(),
                           ebeam));
      SetTreeTitle(Form("IQMD %s events: %s", sim_type.Data(), reac.GetTitle()));
      AddInfo({"Projectile",reac.GetNucleus(KV2Body::projectile)->GetSymbol()});
      AddInfo({"Target",reac.GetNucleus(KV2Body::target)->GetSymbol()});
      AddInfo({"E/A", ebeam});
 
      status = ReadLineAndCheck(3," "); //  0.000 12.000     bmin,bmax
      if(status == ReadStatus::OK)
      {
         AddInfo({"bmin", GetDoubleReadPar(0)});
         AddInfo({"bmax", GetDoubleReadPar(1)});
 
         status = ReadLineAndCheck(4," ");//     1.000  2000.000   1599990     dt,tmax,nevent
         if(status == ReadStatus::OK)
         {
            AddInfo({"dt",GetDoubleReadPar(0)});
            AddInfo({"tmax",GetDoubleReadPar(1)});
            AddInfo({"nevent",GetIntReadPar(2)});
 
            // read 'header' line   "  b(fm)        t    Z    N      px(MeV/n)      py(MeV/n)      pz(MeV/n) event_id"
            return ReadLine()==ReadStatus::OK;
         }
      }
   }
   return false;
}
 
 
 
 
Bool_t KVSimReader_IQMD::ReadEvent(KVSimEvent* evt)
{
   // b(fm)        t    Z    N      px(MeV/n)      py(MeV/n)      pz(MeV/n) event_id
   // 9.518  554.000   69  116     -0.6173435    -10.5668182    -29.6082726 000001
 
   auto status = (reuse_last_line ? ReuseLineAndCheck(8, " ") : ReadLineAndCheck(8," "));
   if(status != ReadStatus::OK)
      return false;
 
   reuse_last_line = false;
   auto event_id = GetIntReadPar(7);
   evt->SetNumber(event_id);
   evt->GetParameters()->SetValue("ImpactParameter", GetDoubleReadPar(0));
   evt->GetParameters()->SetValue("Time", GetDoubleReadPar(1));
   do {
      ReadNucleus(evt->AddParticle());
      status = ReadLineAndCheck(8, " ");
   } while( status==ReadStatus::OK && event_id==GetIntReadPar(7));
 
   if(status == ReadStatus::OK){
      reuse_last_line = true;
   }
 
   return true;
}
 
 
 
 
Bool_t KVSimReader_IQMD::ReadNucleus(KVSimNucleus* nuc)
{
   // Don't actually read anything, just use already read-in values
   // b(fm)        t    Z    N      px(MeV/n)      py(MeV/n)      pz(MeV/n) event_id
   // 9.518  554.000   69  116     -0.6173435    -10.5668182    -29.6082726 000001
 
   nuc->SetZ(GetIntReadPar(2));
   nuc->SetN(GetIntReadPar(3));
   nuc->SetMomentum(nuc->GetA()*GetDoubleReadPar(4),nuc->GetA()*GetDoubleReadPar(5),nuc->GetA()*GetDoubleReadPar(6));
   return true;
}
 
 
ClassImp(KVSimReader_IQMD)