///////////////////////////////////////////////////////////////////////////////
// Hijing2u reads HIJING events and input parameters from the Hijing output
// files, converts them to the UniGen format,
// and saves them to a root file.
//
// James Hostetter & Christoph Baumann (cbaumann@ikf.uni-frankfurt.de)
// March 2010
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include <iomanip>
#include <map>
#include <TROOT.h>
#include <TApplication.h>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TBranch.h>
#include "TMath.h"
#include "UEvent.h"
#include "UParticle.h"
#include "URun.h"
using namespace std;

TROOT root("ROOT", "HIJING to UniGen conversion");

TFile *fi;
TTree *tr;
UEvent *ev;

/*****************************************************************************/
void bomb(char *myst) {
  cerr<<"Error: "<<myst<<", bombing"<<endl;
  exit(-1);
}
/*****************************************************************************/
int trapco(int ityp, int iso3, int ichg) {
  // translate HIJING pid code to pdg code
  static map<Int_t, Int_t> ampt_pdg;
  static int firstcall=1;
  int ampt;
  int pdg;

  if (firstcall) {
    ifstream in;
    const TString unigen = TString(getenv("UNIGEN"));
    //Use UrQMD PDG tables
    in.open((unigen + "/input/urqmd_pdg.dat").Data());
    if (in.fail()) bomb("cannot open urqmd_pdg.dat");
    while(1) {
      if(in.eof()) break;
      in >> ampt >> pdg;
      ampt_pdg[ampt] = pdg;
    }
    in.close();
    firstcall = 0;
  }

  int id;
  if(ityp >= 0) id = 1000 * (ichg+2) + ityp;
  else          id = -1000 * (ichg+2) + ityp;
  if(ampt_pdg.find(id) != ampt_pdg.end()) {
    return ampt_pdg[id];
  }
  return 0;
}
/*****************************************************************************/
int main(int argc, char *argv[]) {
  ifstream in;
  char *inpfile;
  char *outfile;
  char c;
  int nevents;
  string dust;

  URun *ru = 0;
  string version, comment;
  int aproj, zproj, atarg, ztarg, nr, partNum;
  double b, bmin, bmax, sqrts, time, dtime;

  if (argc != 3) {
    //cout << "usage:   " << argv[0] << " HIJING-Output-File outfile nevents\n";
    //cout << "example:   " << argv[0] << " hijing.dat simdat.root 10\n";
    cout << "usage:   " << argv[0] << " HIJING-File ROOT-Outputfile \n";
    cout << "example:   " << argv[0] << " hijing.dat simdat.root\n";
    exit(0);
  }

  inpfile = argv[1];
  outfile = argv[2];
//  nevents = atoi(argv[3]);
  TApplication theApp("App", &argc, argv);

  in.open(inpfile);
  if (in.fail()) bomb("cannot open input file");
  fi = TFile::Open(outfile, "RECREATE", "HIJING");
  tr = new TTree("events","HIJING tree");
  ev = new UEvent();
  tr->Branch("event", "UEvent", &ev);
  
  for (int i=1;i<=52;i++)
  in.ignore(777,'\n');
  //1st Version: Energy read from hijing library output, now: SKIP LINE, for loop does it anyhow
  //  in>>dust>>dust>>sqrts;
  for (int i=1;i<=4;i++)
  in.ignore(777,'\n');
  
  // event loop
  cout << "get global information: " << endl;
  in >> nevents >> dust >> dust >> aproj >> atarg >> zproj >> ztarg >> dust >> bmin >> bmax;
  cout << "Nevents: " << nevents << " A_P: " << aproj << " A_T: " << atarg << " Z_P: " << zproj << " Z_T: " << ztarg << endl;

  //'dust' here holds the Reaction FRAME (CMS/LAB), might be put in proper variable in future
  in >> sqrts >> dust;
  cout<< sqrts<< endl;
  in.ignore(777,'\n');
  const int bunch = 10;
  int events_processed=0;
  for (int n=0; n<nevents; n++) {
    if ((n%bunch)==0) cout << "event "  << setw(5) << n << endl;
    string line;
    char pee;
    pee=in.peek();
    if (pee==EOF) break; 
    //Get all relevant event data from hijing.dat
    in >>dust>> nr >> partNum >> b;

    ev->Clear();
	ev->SetB(b);
    ev->SetEventNr(nr);
    ev->SetPhi(0);
    ev->SetNes((int) (time/dtime));
    events_processed++;

    int step_nr=0;
    double step_time;
    for (int i=0;  i<partNum; i++) {
        double t,x,y,z,e,px,py,pz;
        double weight;
        int ichg,status, parent, parent_decay, mate,ID;
        int decay, child[2];
        in >> ID >> ichg >> weight >> px >> py >> pz;
        if (in.fail()) bomb("while reading tracks");
        status=parent=parent_decay=decay=child[0]=child[1]=0;
        x=y=z=t=0.0;
	e=sqrt(px*px+py*py+pz*pz);
	//weight contains restmass of the particle
        ev->AddParticle(i+1,ID, status, ichg,
                        parent_decay,mate-1, decay,child,
                        px, py, pz, e, x, y, z, t, weight);
      }
      do in.get(c); while (c!='\n');
      ev->SetStepNr(step_nr++);
      ev->SetStepT(step_time);
      /*nout += */tr->Fill();
  }
  cout << events_processed << " events processed\n";
  in.close();

  // create the run object

  string generator = "HIJING";
  generator.append(version);
  double m = 0.938271998;
  double ecm = sqrts/2; // energy per nucleon in cm
  double pcm = sqrt(ecm*ecm-m*m); // momentum per nucleon in cm
  double beta = pcm / ecm;
  double gamma = 1.0/sqrt(1-beta*beta); 
  double pproj = gamma*(+pcm-beta*ecm);
  double ptarg = gamma*(-pcm-beta*ecm);
  ru = new URun(generator.data(), comment.data(), 
      aproj, zproj, pproj, 
      atarg, ztarg, ptarg, 
      bmin, bmax, -1, 0, 0, 0, events_processed);
  ru->Write();
  fi->Write();
  fi->Close();
  return 0;
}
/*****************************************************************************/