Fixed occasional division by zero in epos-tim
[u/mrichter/AliRoot.git] / EPOS / TEpos.cxx
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>
19 #include <TDatabasePDG.h>
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         char *epoEnv = getenv("EPO");
159         if (epoEnv) {
160                 strcpy(epo, epoEnv);
161         } else {
162                 strcpy(epo, getenv("ALICE_ROOT"));
163         }
164         strcat(epo, "/EPOS/epos167");
165
166         file << "fname pathnx " << epo << "/" << endl;
167         file << "fname histo none" << endl;
168         file << "fname copy none" << endl;
169         file << "fname log none" << endl;
170         file << "fname check none" << endl;
171         file << "fname data /tmp/epos.out" << endl;
172         file << "fname initl " << epo << "/epos.initl" << endl;
173         file << "fname inidi " << epo << "/epos.inidi" << endl;
174         file << "fname inidr " << epo << "/epos.inidr" << endl;
175         file << "fname iniev " << epo << "/epos.iniev" << endl;
176         file << "fname inirj " << epo << "/epos.inirj" << endl;
177         file << "fname inics " << epo << "/epos.inics" << endl;
178         file << "fname inigrv " << epo << "/epos.inigrv" << endl;
179         file << "fqgsjet dat " << epo << "/qgsjet/qgsjet.dat" << endl;
180         file << "fqgsjet ncs " << epo << "/qgsjet/qgsjet.ncs" << endl;
181         file << "fqgsjetII dat " << epo << "/qgsjetII/qgsdat-II-03" << endl;
182         file << "fqgsjetII ncs " << epo << "/qgsjetII/sectnu-II-03" << endl;
183         file << "nodecay 120" << endl;
184         file << "nodecay -120" << endl;
185         file << "nodecay 130" << endl;
186         file << "nodecay -130" << endl;
187         file << "nodecay -20" << endl;
188         file << "nodecay 20" << endl;
189         file << "nodecay 14" << endl;
190         file << "nodecay -14" << endl;
191         file << "set ndecay 1111110" << endl;
192         file << "echo on" << endl;
193
194 //      .optns file
195         int precision = file.precision();
196         file.precision(15);
197         file << "set seedj " << (gRandom->Rndm() * 1e14) << endl;
198         file.precision(precision);
199         file << "application hadron" << endl;
200         file << "set laproj " << fLaproj << endl;
201         file << "set maproj " << fMaproj << endl;
202         file << "set latarg " << fLatarg << endl;
203         file << "set matarg " << fMatarg << endl;
204         file << "set bminim " << fBminim << endl;
205         file << "set bmaxim " << fBmaxim << endl;
206         file << "set phimin " << fPhimin << endl;
207         file << "set phimax " << fPhimax << endl;
208         file << "set ecms " << fEcms << endl;
209
210         for(unsigned int i = 0; i < fNoDecays.size(); ++i) {
211                 file << "nodecay " << fNoDecays[i] << endl;
212         }
213
214         file << "switch splitting " << (fSplitting ? "on" : "off") << endl;
215
216         file << "frame nucleon-nucleon" << endl;
217
218         for(unsigned int i = 0; i < fExtraInputLines.size(); ++i) {
219                 file << fExtraInputLines[i] << endl;
220         }
221
222     file << "output epos" << endl;
223     file << "record event nevt nptl b endrecord" << endl;
224     file << "record particle i id fa mo c1 c2 st endrecord" << endl;
225
226         file << "input " << epo << "/epos.param" << endl;
227         file << "runprogram" << endl;
228         file.close();
229 }
230
231 Float_t TEpos::GetBimevt() { return cevt.bimevt; }
232 Float_t TEpos::GetPhievt() { return cevt.phievt; }
233 Int_t TEpos::GetKolevt() { return cevt.kolevt; }
234 Int_t TEpos::GetKoievt() { return cevt.koievt; }
235 Float_t TEpos::GetPmxevt() { return cevt.pmxevt; }
236 Float_t TEpos::GetEgyevt() { return cevt.egyevt; }
237 Int_t TEpos::GetNpjevt() { return cevt.npjevt; }
238 Int_t TEpos::GetNtgevt() { return cevt.ntgevt; }
239 Int_t TEpos::GetNpnevt() { return cevt.npnevt; }
240 Int_t TEpos::GetNppevt() { return cevt.nppevt; }
241 Int_t TEpos::GetNtnevt() { return cevt.ntnevt; }
242 Int_t TEpos::GetNtpevt() { return cevt.ntpevt; }
243 Int_t TEpos::GetJpnevt() { return cevt.jpnevt; }
244 Int_t TEpos::GetJppevt() { return cevt.jppevt; }
245 Int_t TEpos::GetJtnevt() { return cevt.jtnevt; }
246 Int_t TEpos::GetJtpevt() { return cevt.jtpevt; }
247 Float_t TEpos::GetXbjevt() { return cevt.xbjevt; }
248 Float_t TEpos::GetQsqevt() { return cevt.qsqevt; }
249 Int_t TEpos::GetNglevt() { return cevt.nglevt; }
250 Float_t TEpos::GetZppevt() { return cevt.zppevt; }
251 Float_t TEpos::GetZptevt() { return cevt.zptevt; }
252