Просмотр исходного кода

Update urqmd2u to meet the new UDst version.

Nikita Ermakov лет назад: 5
Родитель
Сommit
d3a4310d55
1 измененных файлов с 176 добавлено и 167 удалено
  1. 176 167
      urqmd2u.cpp

+ 176 - 167
urqmd2u.cpp

@@ -24,196 +24,205 @@
 #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;
-UEvent *ev;
+TClonesArray *arrays[UArrays::NAllUArrays];
 
 //_________________
 void bomb(const char *myst) {
-  std::cerr << "Error: " << myst << ", bombing" << std::endl;
-  exit(-1);
+	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);
+	// 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;
-  std::string dust;
-
-  URun *ru = 0;
-  std::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) {
-    std::cout << "usage:   " << argv[0] << " inpfile outfile nevents\n";
-    std::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;
-
-  tr = new TTree("events","UrQMD tree", split);
-  tr->SetAutoSave(1000000);
-  ev = new UEvent();
-  tr->Branch("event", "UEvent", &ev, bufsize, split);
-
-  // event loop
-
-  const int bunch = 10;
-
-  int events_processed=0;
-  for (int n=0; n<nevents; n++) {
-    if ((n%bunch)==0) std::cout << "event "  << setw(5) << n << std::endl;
-    std::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
-
-    std::cout << "version: " << version << " sqrts: " << sqrts << " dtime: " << dtime << std::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'); 
-
-    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;
-      
-      std::cout << "Number of particles in event: " << mult << std::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++) {
-	std::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;
+	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);
 	}
-	std::cout << "px: " << px << " py: " << py << " pz: " << pz << std::endl;
-	if (in.fail()) bomb("while reading tracks");
-	status=parent=parent_decay=decay=child[0]=child[1]=0;
-	std::cout << "I am here 1" << std::endl;
-	ev->addParticle(i, trapco(ityp, ichg), status, parent,
-			parent_decay, mate-1, decay, child,
-			px, py, pz, e, x, y, z, t, weight);
-	std::cout << "I am here 2" << std::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();
-  std::cout << events_processed << " events processed\n";
-
-  // create the run object
-
-  std::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;
+
+	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;
 }