unigen/urqmd2u.cpp

///////////////////////////////////////////////////////////////////////////////
// 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
///////////////////////////////////////////////////////////////////////////////

// C++ headers
#include <iostream>
#include <fstream>
#include <iomanip>
#include <map>

// ROOT headers
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TBranch.h"
#include "TMath.h"
#include "TClonesArray.h"

// UniGen headers
#include "URun.h"
#include "UEvent.h"
#include "UParticle.h"
#include "UPIDConverter.h"
#include "UArrays.h"

using namespace std;

TFile *fi;
TTree *tr;
TClonesArray *arrays[UArrays::NAllUArrays];

//_________________
void bomb(const char *myst) {
	cerr << "Error: " << myst << ", bombing" << endl;
	exit(-1);
}

//_________________
int trapco(int ityp, int ichg) {
	// translate UrQMD pid code to pdg code

	/* UrQMD PIDs are in fact composite - a particle is fully defined by the
	   type specifier (ityp), the charge (ichg) and in case of baryons, the
	   third component isospin (iso3; ignored by us at present). For
	   simplicity, our conversion tables collapse these into a single number
	   as follows:
	   - propagate the sign of ityp (particle-antiparticle distinction for
	   baryons, strangeness-antistrangeness distinction for mesons) to that
	   of the stored value;
	   - shift the ichg range from -2..2 (UrQMD does not support nuclear
	   fragments other than protons and neutrons so all particles it
	   produces fall in this range) to 0..4 to make sure it doesn't
	   interfere with the above;
	   - multiply shifted charge by +/-1000 and add it to type. The latter
	   is guaranteed to be smaller than 1000 (baryon types are one- or
	   two-digit, meson types three-digit) so that way no ambiguities
	   occur. */
	int id;
	if (ityp >= 0) id = 1000 * (ichg + 2) + ityp;
	else
		id = -1000 * (ichg + 2) + ityp;

	/* This will return 0 for unknown input values. */
	return UPIDConverter::Instance()->GetPDGCode(id, UPIDConverter::eUrQMD);
}

//_________________
int main(int argc, char *argv[]) {

	ifstream in;
	char *inpfile;
	char *outfile;
	char c;
	int nevents;
	string dust;

	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]);

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

	fi = TFile::Open(outfile, "RECREATE", "UrQMD");
	fi->SetCompressionLevel(9);
	int bufsize = 65536;
	int split = 99;

	// Create and set out TTree
	tr = new TTree("events","UrQMD tree", split);
	tr->SetAutoSave(1000000);
	for (unsigned int i = 0; i < UArrays::NAllUArrays; i++) {
		// Create arrayss
		arrays[i] = new TClonesArray(UArrays::uArrayTypes[i],
									 UArrays::uArraySizes[i]);
		// Set branch
		tr->Branch(UArrays::uArrayNames[i], &arrays[i], bufsize, split);
	}

	// 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

		cout << "version: " << version << " sqrts: " << sqrts << " dtime: " << dtime << endl;
		comment.clear();
		// read 4 lines of options and 6 lines of params
		for (int i=0; i<10; i++) {
			getline(in,line); 
			comment.append(line);
			comment.append("\n");
		}
		in.ignore(777,'\n'); 

		// Let the new UEvent appear!
		UEvent *ev = new ((*(arrays[UArrays::Event]))[arrays[UArrays::Event]->GetEntries()]) UEvent();

		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;

			cout << "Number of particles in event: " << mult << endl;

			in.ignore(777,'\n'); // ignore the rest of the line
			getline(in,line);
			ev->setComment(line.data());

			for (int i=0;  i<mult; i++) {
				cout << "Working on particle i: " << i;
				double t,x,y,z,e,px,py,pz;
				double weight = 1.0;
				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;
				}
				cout << "px: " << px << " py: " << py << " pz: " << pz << endl;
				if (in.fail()) bomb("while reading tracks");
				status=parent=parent_decay=decay=child[0]=child[1]=0;
				cout << "I am here 1" << endl;
				new ((*(arrays[UArrays::Particle]))[arrays[UArrays::Particle]->GetEntries()]) UParticle(i, trapco(ityp, ichg), status, parent,
																										parent_decay, mate-1, decay, child,
																										px, py, pz, e, x, y, z, t, weight);
				cout << "I am here 2" << endl;
			} // for (int i=0;  i<mult; i++)

			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);
	new ((*(arrays[UArrays::Run]))[arrays[UArrays::Run]->GetEntries()]) URun(generator.data(), comment.data(),
																			 aproj, zproj, pproj, 
																			 atarg, ztarg, ptarg, 
																			 bmin, bmax, -1, 0, 0, sigma, events_processed);
	nout += tr->Fill();
	fi->Write();
	fi->Close();
	cout << "Total bytes were written: " << nout << endl;
	return 0;
}