///////////////////////////////////////////////////////////////////////////////
// urqmd2u reads UrQMD events from the ftn13 or ftn14 ascii files and
// converts them to the UniGen format and saves on a root file.
//
// ftn14 contains snapshots at given times (event steps). The event steps 
// are converted to separate events. 
//
// ftn13 contains the final snapshot and the freeze-out coordinates. 
// The freeze-out coordinates are used. The final snapshot coordinates 
// are discarded.
//
// D.Miskowiec, February 2006
///////////////////////////////////////////////////////////////////////////////
#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 "URun.h"
#include "UEvent.h"
#include "UParticle.h"
using namespace std;

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

//TROOT root("ROOT", "UrQMD to UniGen conversion");

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

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

  int id;
  if(ityp >= 0) id = 1000 * (ichg+2) + ityp;
  else          id = -1000 * (ichg+2) + ityp;
  if(urqmd_pdg.find(id) != urqmd_pdg.end()) {
      return urqmd_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 filetype, eos, aproj, zproj, atarg, ztarg, nr;
  double beta, b, bmin, bmax, sigma, elab, plab, sqrts, time, dtime;

  if (argc != 4) {
    cout << "usage:   " << argv[0] << " inpfile outfile nevents\n";
    cout << "example: " << argv[0] << " ftn14 ftn14.root 10\n";
    exit(0);
  }

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

  int nout=0;
  in.open(inpfile);
  if (in.fail()) bomb("cannot open input file");

  fi = TFile::Open(outfile, "RECREATE", "UrQMD");

  tr = new TTree("events","UrQMD tree");
  ev = new UEvent();
  tr->Branch("event", "UEvent", &ev);

  // event loop

  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;

    in >> dust >> dust >> version >> dust >> dust;
    in >> dust >> filetype >> dust >> dust >> dust >> aproj >> zproj;
    in >> dust >> dust >> dust >> atarg >> ztarg; 
    in >> dust >> dust >> dust >> beta >> dust >> dust;
    in >> dust >> b >> bmin >> bmax >> dust >> sigma;
    in >> dust >> eos >> dust >> elab >> dust >> sqrts >> dust >> plab;
    in >> dust >> nr >> dust >> dust >> dust;
    in >> dust >> dust >> time >> dust >> dtime;
    in.ignore(777,'\n'); // ignore the rest of the line

    comment.clear();
    // read 3 lines of options and 4 lines of params
    for (int i=0; i<7; i++) {
      getline(in,line); 
      comment.append(line);
      comment.append("\n");
    }
    in.ignore(777,'\n'); 

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

    int step_nr=0;
    char pee;
    while (1) { // loop over time slices
      int mult;
      double step_time;
      pee=in.peek();
      if (pee=='U') break; 
      if (pee==EOF) break; 
      in >> mult >> step_time; 
      in.ignore(777,'\n'); // ignore the rest of the line
      getline(in,line);
      ev->SetComment(line.data());
      for (int i=0;  i<mult; i++) {
	double t,x,y,z,e,px,py,pz;
	double weight;
	int ityp, iso3, ichg, status, parent, parent_decay, mate;
	int decay, child[2];
	in >> t >> x >> y >> z;
	in >> e >> px >> py >> pz >> dust;
	in >> ityp >> iso3 >> ichg >> mate >> dust >> dust;
	if (filetype==13) { // freeze-out coordinates
	  in >> t >> x >> y >> z;
	  in >> e >> px >> py >> pz;
	}
	if (in.fail()) bomb("while reading tracks");
	status=parent=parent_decay=decay=child[0]=child[1]=0;
	weight = 1.0;
	ev->AddParticle(i, trapco(ityp,/*iso3,*/ichg), status, parent,
			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();
    }
    if (pee==EOF) break; 
  }
  in.close();
  cout << events_processed << " events processed\n";

  // create the run object

  string generator = "UrQMD";
  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 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, sigma, events_processed);
  ru->Write();
  fi->Write();
  fi->Close();
  return nout;
  return 0;
}
/*****************************************************************************/