[u/mrichter/AliRoot.git] / EPOS / TEpos.cxx

/*
 *###################################################################
 *#        EPOS 1.67     K. WERNER, T. PIEROG, S. PORTEBOEUF.       #
 *#                      Contact: werner@subatech.in2p3.fr          #
 *###################################################################
 *
 * TEpos.cxx
 * 
 * Wraper class for interfacing EPOS model, derived from ROOT's TGenerator.
 * It generates temporary input file for the model, providing user with
 * ability to add his/hers own lines to the input.
 * Output is read directly from common blocks.
 *
 *      Author: Piotr Ostrowski, postrow@if.pw.edu.pl
 */


#include <TClonesArray.h>
#include <TDatabasePDG.h>
#include <TObjArray.h>
#include <TParticle.h>
#include <TROOT.h>
#include <TRandom.h>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "eposproc.h"
#include "EPOScommon.h"
#include "TEpos.h"

using namespace std;

ClassImp(TEpos)

TEpos::TEpos()  : TGenerator("EPOS", "Epos event generator"),
		fLaproj(0),
		fMaproj(0),
		fLatarg(0),
		fMatarg(0),
		fBminim(0.0),
		fBmaxim(10000.0),
		fPhimin(0.0),
		fPhimax(2*3.1415927),
		fEcms(-1),
		fSplitting(kFALSE),
		fNoDecays(),
		fExtraInputLines()
{
}

TEpos::~TEpos() {}

void TEpos::Initialize() {
	Int_t nopeno = 0;
	GenerateInputFile();
	aaset(0);
	atitle();
	xiniall();

	const char *inputFileName = GetInputFileName();
	Int_t nameLength = strlen(inputFileName);
	setinp(inputFileName, nameLength, nameLength);
	aread();
	while(copen.nopen == -1) {
		copen.nopen=nopeno;
	    prnt1.iecho=1;
		xiniall();
	    aread();
	}
	Int_t ishini;
	utpri("aamain",prnt1.ish,ishini,4,6);
	if(appli.model != 1)
		IniModel(appli.model);
	ebin.nrebin = 1;
	ainit();
}

void TEpos::GenerateEvent() {
	cseed.seedj = gRandom->Rndm() * 1e10;
	aseed(2);
	Int_t n = 1;
	evgen(n);
}

Int_t TEpos::ImportParticles(TClonesArray *particles, Option_t *) {
	particles->Clear();
	if (!cevt.nevt) return 0;
	Int_t numpart = cptl.nptl;
	printf("%d particles generated\n", numpart);
	for (Int_t iPart=0; iPart<numpart; iPart++) {
		Int_t tFather = cptl.iorptl[iPart] - 1;
		tFather = tFather < -1 ? -1 : tFather;
		if (tFather> -1) {
			TParticle *mother = (TParticle*) (particles->UncheckedAt(tFather));
      			mother->SetLastDaughter(iPart);
      			if (mother->GetFirstDaughter()==-1)
 				mother->SetFirstDaughter(iPart);
    		}
		TDatabasePDG *pdgDb = TDatabasePDG::Instance();
		Int_t status = cptl.istptl[iPart] + 1;
		Int_t pdg = pdgDb->ConvertIsajetToPdg(cptl.idptl[iPart]);
		if (pdg == 0) {
			printf("TEpos: Warning, unknown particle, index: %d, ISAJET: %d\n",iPart,cptl.idptl[iPart]);
		}
    		new((*particles)[iPart]) TParticle(pdg, status,
				 tFather, -1, -1, -1,
				 cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
				 cptl.xorptl[iPart][0]*1.e-12, cptl.xorptl[iPart][1]*1.e-12, cptl.xorptl[iPart][2]*1.e-12, cptl.xorptl[iPart][3]*1e-12);
    		(*particles)[iPart]->SetUniqueID(iPart);
  	}

  	return numpart;
}

TObjArray*  TEpos::ImportParticles(Option_t *) {
	fParticles->Clear();
	if (!cevt.nevt) return NULL;
	Int_t numpart = cptl.nptl;
	printf("%d particles generated\n", numpart);
	for (Int_t iPart=0; iPart<numpart; iPart++) {
		Int_t tFather = cptl.iorptl[iPart] - 1;
		tFather = tFather < -1 ? -1 : tFather;
		if (tFather> -1) {
			TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
      			mother->SetLastDaughter(iPart);
      			if (mother->GetFirstDaughter()==-1)
 				mother->SetFirstDaughter(iPart);
    		}
		TDatabasePDG *pdgDb = TDatabasePDG::Instance();
		Int_t status = cptl.istptl[iPart] + 1;
		Int_t pdg = pdgDb->ConvertIsajetToPdg(cptl.idptl[iPart]);
		if (pdg == 0) {
			printf("TEpos: Warning, unknown particle, index: %d, ISAJET: %d\n",iPart,cptl.idptl[iPart]);
		}
    		TParticle* p = new TParticle(pdg, status,
				 tFather, -1, -1, -1,
				 cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
				 cptl.xorptl[iPart][0]*1.e-12, cptl.xorptl[iPart][1]*1.e-12, cptl.xorptl[iPart][2]*1.e-12, cptl.xorptl[iPart][3]*1e-12);
    		p->SetUniqueID(iPart);
    		fParticles->Add(p);
  	}

  	return fParticles;
}

void TEpos::AddNoDecay(Int_t nodecay) {
	fNoDecays.push_back(nodecay);
}

void TEpos::AddExtraInputLine(const char *line) {
	fExtraInputLines.push_back(line);
}

void TEpos::GenerateInputFile() {
	ofstream file(GetInputFileName(), ios_base::out | ios_base::trunc);
	char epo[256];
	strcpy(epo, getenv("ALICE_ROOT"));
	strcat(epo, "/EPOS/epos167");

	file << "fname pathnx " << epo << "/" << endl;
	file << "fname histo none" << endl;
	file << "fname copy none" << endl;
	file << "fname log none" << endl;
	file << "fname check none" << endl;
	file << "fname data /tmp/epos.out" << endl;
	file << "fname initl " << epo << "/epos.initl" << endl;
	file << "fname inidi " << epo << "/epos.inidi" << endl;
	file << "fname inidr " << epo << "/epos.inidr" << endl;
	file << "fname iniev " << epo << "/epos.iniev" << endl;
	file << "fname inirj " << epo << "/epos.inirj" << endl;
	file << "fname inics " << epo << "/epos.inics" << endl;
	file << "fname inigrv " << epo << "/epos.inigrv" << endl;
	file << "fqgsjet dat " << epo << "/qgsjet/qgsjet.dat" << endl;
	file << "fqgsjet ncs " << epo << "/qgsjet/qgsjet.ncs" << endl;
	file << "fqgsjetII dat " << epo << "/qgsjetII/qgsdat-II-03" << endl;
	file << "fqgsjetII ncs " << epo << "/qgsjetII/sectnu-II-03" << endl;
	file << "nodecay 120" << endl;
	file << "nodecay -120" << endl;
	file << "nodecay 130" << endl;
	file << "nodecay -130" << endl;
	file << "nodecay -20" << endl;
	file << "nodecay 20" << endl;
	file << "nodecay 14" << endl;
	file << "nodecay -14" << endl;
	file << "set ndecay 1111110" << endl;
	file << "echo on" << endl;

//	.optns file
	file << "application hadron" << endl;
	file << "set laproj " << fLaproj << endl;
	file << "set maproj " << fMaproj << endl;
	file << "set latarg " << fLatarg << endl;
	file << "set matarg " << fMatarg << endl;
	file << "set bminim " << fBminim << endl;
	file << "set bmaxim " << fBmaxim << endl;
	file << "set phimin " << fPhimin << endl;
	file << "set phimax " << fPhimax << endl;
	file << "set ecms " << fEcms << endl;

	for(unsigned int i = 0; i < fNoDecays.size(); ++i) {
		file << "nodecay " << fNoDecays[i] << endl;
	}

	file << "switch splitting " << (fSplitting ? "on" : "off") << endl;

	file << "frame nucleon-nucleon" << endl;

	for(unsigned int i = 0; i < fExtraInputLines.size(); ++i) {
		file << fExtraInputLines[i] << endl;
	}

    file << "output epos" << endl;
    file << "record event nevt nptl b endrecord" << endl;
    file << "record particle i id fa mo c1 c2 st endrecord" << endl;

	file << "input " << epo << "/epos.param" << endl;
	file << "runprogram" << endl;
	file.close();
}

Float_t TEpos::GetBimevt() { return cevt.bimevt; }
Float_t TEpos::GetPhievt() { return cevt.phievt; }
Int_t TEpos::GetKolevt() { return cevt.kolevt; }
Int_t TEpos::GetKoievt() { return cevt.koievt; }
Float_t TEpos::GetPmxevt() { return cevt.pmxevt; }
Float_t TEpos::GetEgyevt() { return cevt.egyevt; }
Int_t TEpos::GetNpjevt() { return cevt.npjevt; }
Int_t TEpos::GetNtgevt() { return cevt.ntgevt; }
Int_t TEpos::GetNpnevt() { return cevt.npnevt; }
Int_t TEpos::GetNppevt() { return cevt.nppevt; }
Int_t TEpos::GetNtnevt() { return cevt.ntnevt; }
Int_t TEpos::GetNtpevt() { return cevt.ntpevt; }
Int_t TEpos::GetJpnevt() { return cevt.jpnevt; }
Int_t TEpos::GetJppevt() { return cevt.jppevt; }
Int_t TEpos::GetJtnevt() { return cevt.jtnevt; }
Int_t TEpos::GetJtpevt() { return cevt.jtpevt; }
Float_t TEpos::GetXbjevt() { return cevt.xbjevt; }
Float_t TEpos::GetQsqevt() { return cevt.qsqevt; }
Int_t TEpos::GetNglevt() { return cevt.nglevt; }
Float_t TEpos::GetZppevt() { return cevt.zppevt; }
Float_t TEpos::GetZptevt() { return cevt.zptevt; }
Commit	Line	Data
9ef1c2d9	1	/*
	2	*###################################################################
	3	*# EPOS 1.67 K. WERNER, T. PIEROG, S. PORTEBOEUF. #
	4	*# Contact: werner@subatech.in2p3.fr #
	5	*###################################################################
	6	*
	7	* TEpos.cxx
	8	*
	9	* Wraper class for interfacing EPOS model, derived from ROOT's TGenerator.
	10	* It generates temporary input file for the model, providing user with
	11	* ability to add his/hers own lines to the input.
	12	* Output is read directly from common blocks.
	13	*
	14	* Author: Piotr Ostrowski, postrow@if.pw.edu.pl
	15	*/
	16
	17
	18	#include <TClonesArray.h>
6a0522a4	19	#include <TDatabasePDG.h>
9ef1c2d9	20	#include <TObjArray.h>
	21	#include <TParticle.h>
	22	#include <TROOT.h>
	23	#include <TRandom.h>
	24	#include <iostream>
	25	#include <fstream>
	26	#include <string>
	27	#include <vector>
	28	#include "eposproc.h"
	29	#include "EPOScommon.h"
	30	#include "TEpos.h"
	31
	32	using namespace std;
	33
	34	ClassImp(TEpos)
	35
	36	TEpos::TEpos() : TGenerator("EPOS", "Epos event generator"),
	37	fLaproj(0),
	38	fMaproj(0),
	39	fLatarg(0),
	40	fMatarg(0),
	41	fBminim(0.0),
	42	fBmaxim(10000.0),
	43	fPhimin(0.0),
	44	fPhimax(2*3.1415927),
	45	fEcms(-1),
	46	fSplitting(kFALSE),
	47	fNoDecays(),
	48	fExtraInputLines()
	49	{
	50	}
	51
	52	TEpos::~TEpos() {}
	53
	54	void TEpos::Initialize() {
	55	Int_t nopeno = 0;
	56	GenerateInputFile();
	57	aaset(0);
	58	atitle();
	59	xiniall();
	60
	61	const char *inputFileName = GetInputFileName();
	62	Int_t nameLength = strlen(inputFileName);
	63	setinp(inputFileName, nameLength, nameLength);
	64	aread();
	65	while(copen.nopen == -1) {
	66	copen.nopen=nopeno;
	67	prnt1.iecho=1;
	68	xiniall();
	69	aread();
	70	}
	71	Int_t ishini;
	72	utpri("aamain",prnt1.ish,ishini,4,6);
	73	if(appli.model != 1)
	74	IniModel(appli.model);
	75	ebin.nrebin = 1;
	76	ainit();
	77	}
	78
	79	void TEpos::GenerateEvent() {
	80	cseed.seedj = gRandom->Rndm() * 1e10;
	81	aseed(2);
	82	Int_t n = 1;
	83	evgen(n);
84	}
85
86	Int_t TEpos::ImportParticles(TClonesArray particles, Option_t ) {
87	particles->Clear();
88	if (!cevt.nevt) return 0;
89	Int_t numpart = cptl.nptl;
90	printf("%d particles generated\n", numpart);
91	for (Int_t iPart=0; iPart<numpart; iPart++) {
92	Int_t tFather = cptl.iorptl[iPart] - 1;
93	tFather = tFather < -1 ? -1 : tFather;
94	if (tFather> -1) {
95	TParticle mother = (TParticle) (particles->UncheckedAt(tFather));
96	mother->SetLastDaughter(iPart);
97	if (mother->GetFirstDaughter()==-1)
98	mother->SetFirstDaughter(iPart);
99	}
100	TDatabasePDG *pdgDb = TDatabasePDG::Instance();
101	Int_t status = cptl.istptl[iPart] + 1;
102	Int_t pdg = pdgDb->ConvertIsajetToPdg(cptl.idptl[iPart]);
103	if (pdg == 0) {
104	printf("TEpos: Warning, unknown particle, index: %d, ISAJET: %d\n",iPart,cptl.idptl[iPart]);
105	}
106	new((*particles)[iPart]) TParticle(pdg, status,
107	tFather, -1, -1, -1,
108	cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
109	cptl.xorptl[iPart][0]1.e-12, cptl.xorptl[iPart][1]1.e-12, cptl.xorptl[iPart][2]1.e-12, cptl.xorptl[iPart][3]1e-12);
110	(*particles)[iPart]->SetUniqueID(iPart);
111	}
112
113	return numpart;
114	}
115
116	TObjArray* TEpos::ImportParticles(Option_t *) {
117	fParticles->Clear();
118	if (!cevt.nevt) return NULL;
119	Int_t numpart = cptl.nptl;
120	printf("%d particles generated\n", numpart);
121	for (Int_t iPart=0; iPart<numpart; iPart++) {
122	Int_t tFather = cptl.iorptl[iPart] - 1;
123	tFather = tFather < -1 ? -1 : tFather;
124	if (tFather> -1) {
125	TParticle mother = (TParticle) (fParticles->UncheckedAt(tFather));
126	mother->SetLastDaughter(iPart);
127	if (mother->GetFirstDaughter()==-1)
128	mother->SetFirstDaughter(iPart);
129	}
130	TDatabasePDG *pdgDb = TDatabasePDG::Instance();
131	Int_t status = cptl.istptl[iPart] + 1;
132	Int_t pdg = pdgDb->ConvertIsajetToPdg(cptl.idptl[iPart]);
133	if (pdg == 0) {
134	printf("TEpos: Warning, unknown particle, index: %d, ISAJET: %d\n",iPart,cptl.idptl[iPart]);
135	}
136	TParticle* p = new TParticle(pdg, status,
137	tFather, -1, -1, -1,
138	cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
139	cptl.xorptl[iPart][0]1.e-12, cptl.xorptl[iPart][1]1.e-12, cptl.xorptl[iPart][2]1.e-12, cptl.xorptl[iPart][3]1e-12);
140	p->SetUniqueID(iPart);
141	fParticles->Add(p);
142	}
143
144	return fParticles;
145	}
146
147	void TEpos::AddNoDecay(Int_t nodecay) {
148	fNoDecays.push_back(nodecay);
149	}
150
151	void TEpos::AddExtraInputLine(const char *line) {
152	fExtraInputLines.push_back(line);
153	}
154
155	void TEpos::GenerateInputFile() {
156	ofstream file(GetInputFileName(), ios_base::out \| ios_base::trunc);
157	char epo[256];
158	strcpy(epo, getenv("ALICE_ROOT"));
159	strcat(epo, "/EPOS/epos167");
160
161	file << "fname pathnx " << epo << "/" << endl;
162	file << "fname histo none" << endl;
163	file << "fname copy none" << endl;
164	file << "fname log none" << endl;
165	file << "fname check none" << endl;
166	file << "fname data /tmp/epos.out" << endl;
167	file << "fname initl " << epo << "/epos.initl" << endl;
168	file << "fname inidi " << epo << "/epos.inidi" << endl;
169	file << "fname inidr " << epo << "/epos.inidr" << endl;
170	file << "fname iniev " << epo << "/epos.iniev" << endl;
171	file << "fname inirj " << epo << "/epos.inirj" << endl;
172	file << "fname inics " << epo << "/epos.inics" << endl;
173	file << "fname inigrv " << epo << "/epos.inigrv" << endl;
174	file << "fqgsjet dat " << epo << "/qgsjet/qgsjet.dat" << endl;
175	file << "fqgsjet ncs " << epo << "/qgsjet/qgsjet.ncs" << endl;
176	file << "fqgsjetII dat " << epo << "/qgsjetII/qgsdat-II-03" << endl;
177	file << "fqgsjetII ncs " << epo << "/qgsjetII/sectnu-II-03" << endl;
178	file << "nodecay 120" << endl;
179	file << "nodecay -120" << endl;
180	file << "nodecay 130" << endl;
181	file << "nodecay -130" << endl;
182	file << "nodecay -20" << endl;
183	file << "nodecay 20" << endl;
184	file << "nodecay 14" << endl;
185	file << "nodecay -14" << endl;
186	file << "set ndecay 1111110" << endl;
187	file << "echo on" << endl;
188
189	// .optns file
190	file << "application hadron" << endl;
191	file << "set laproj " << fLaproj << endl;
192	file << "set maproj " << fMaproj << endl;
193	file << "set latarg " << fLatarg << endl;
194	file << "set matarg " << fMatarg << endl;
195	file << "set bminim " << fBminim << endl;
196	file << "set bmaxim " << fBmaxim << endl;
197	file << "set phimin " << fPhimin << endl;
198	file << "set phimax " << fPhimax << endl;
199	file << "set ecms " << fEcms << endl;
200
201	for(unsigned int i = 0; i < fNoDecays.size(); ++i) {
202	file << "nodecay " << fNoDecays[i] << endl;
203	}
204
205	file << "switch splitting " << (fSplitting ? "on" : "off") << endl;
206
207	file << "frame nucleon-nucleon" << endl;
208
209	for(unsigned int i = 0; i < fExtraInputLines.size(); ++i) {
210	file << fExtraInputLines[i] << endl;
211	}
212
213	file << "output epos" << endl;
214	file << "record event nevt nptl b endrecord" << endl;
215	file << "record particle i id fa mo c1 c2 st endrecord" << endl;
216
217	file << "input " << epo << "/epos.param" << endl;
218	file << "runprogram" << endl;
219	file.close();
220	}
221
222	Float_t TEpos::GetBimevt() { return cevt.bimevt; }
223	Float_t TEpos::GetPhievt() { return cevt.phievt; }
224	Int_t TEpos::GetKolevt() { return cevt.kolevt; }
225	Int_t TEpos::GetKoievt() { return cevt.koievt; }
226	Float_t TEpos::GetPmxevt() { return cevt.pmxevt; }
227	Float_t TEpos::GetEgyevt() { return cevt.egyevt; }
228	Int_t TEpos::GetNpjevt() { return cevt.npjevt; }
229	Int_t TEpos::GetNtgevt() { return cevt.ntgevt; }
230	Int_t TEpos::GetNpnevt() { return cevt.npnevt; }
231	Int_t TEpos::GetNppevt() { return cevt.nppevt; }
232	Int_t TEpos::GetNtnevt() { return cevt.ntnevt; }
233	Int_t TEpos::GetNtpevt() { return cevt.ntpevt; }
234	Int_t TEpos::GetJpnevt() { return cevt.jpnevt; }
235	Int_t TEpos::GetJppevt() { return cevt.jppevt; }
236	Int_t TEpos::GetJtnevt() { return cevt.jtnevt; }
237	Int_t TEpos::GetJtpevt() { return cevt.jtpevt; }
238	Float_t TEpos::GetXbjevt() { return cevt.xbjevt; }
239	Float_t TEpos::GetQsqevt() { return cevt.qsqevt; }
240	Int_t TEpos::GetNglevt() { return cevt.nglevt; }
241	Float_t TEpos::GetZppevt() { return cevt.zppevt; }
242	Float_t TEpos::GetZptevt() { return cevt.zptevt; }
243