Reading muon trigger scalers with the DA of the muon trigger and transfer
[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 <TMath.h>
25 #include <iostream>
26 #include <fstream>
27 #include <string>
28 #include <vector>
29 #include "eposproc.h"
30 #include "EPOScommon.h"
31 #include "AliGenEposIsajetToPdgConverter.h"
32 #include "TEpos.h"
33
34 using namespace std;
35
36 ClassImp(TEpos)
37
38 TEpos::TEpos()  : TGenerator("EPOS", "Epos event generator"),
39                 fLaproj(0),
40                 fMaproj(0),
41                 fLatarg(0),
42                 fMatarg(0),
43                 fBminim(0.0),
44                 fBmaxim(10000.0),
45                 fPhimin(0.0),
46                 fPhimax(TMath::TwoPi()),
47                 fEcms(-1),
48                 fSplitting(kFALSE),
49                 fNoDecays(),
50                 fExtraInputLines(),
51                 fIdConverter()
52 {
53         fIdConverter=new AliGenEposIsajetToPdgConverter();
54 }
55
56 TEpos::TEpos(const TEpos&)  : TGenerator("EPOS", "Epos event generator"),
57                 fLaproj(0),
58                 fMaproj(0),
59                 fLatarg(0),
60                 fMatarg(0),
61                 fBminim(0.0),
62                 fBmaxim(10000.0),
63                 fPhimin(0.0),
64                 fPhimax(TMath::TwoPi()),
65                 fEcms(-1),
66                 fSplitting(kFALSE),
67                 fNoDecays(),
68                 fExtraInputLines(),
69                 fIdConverter()
70 {
71         fIdConverter = new AliGenEposIsajetToPdgConverter();
72 }
73
74 TEpos::~TEpos() {
75         delete fIdConverter;
76 }
77
78 TEpos& TEpos::operator=(const TEpos& epos) {
79   //operator=
80   if (this != &epos) {
81     if (fIdConverter) {
82         delete fIdConverter;
83         fIdConverter = new AliGenEposIsajetToPdgConverter(*epos.fIdConverter);
84     } else{
85         fIdConverter = 0;
86     }
87     fLaproj = epos.fLaproj;
88     fMaproj = epos.fMaproj;
89     fLatarg = epos.fLatarg;
90     fMatarg = epos.fMatarg;
91     fBminim = epos.fBminim;
92     fBmaxim = epos.fBmaxim;
93     fPhimin = epos.fPhimin;
94     fPhimax = epos.fPhimax;
95     fEcms   = epos.fEcms;
96     fSplitting = epos.fSplitting;
97     for(UInt_t i = 0; i < epos.fNoDecays.size(); ++i) {
98       fNoDecays.push_back(epos.fNoDecays[i]);
99     }
100     for(UInt_t i = 0; i < epos.fExtraInputLines.size(); ++i) {
101       fExtraInputLines.push_back(epos.fExtraInputLines[i]);
102     }
103   }
104   return *this;
105 }
106
107 void TEpos::Initialize() {
108   // Generates input file and prepares EPOS to read from it.
109         Int_t nopeno = 0;
110         GenerateInputFile();
111         aaset(0);
112         atitle();
113         xiniall();
114
115         const char *inputFileName = GetInputFileName();
116         Int_t nameLength = strlen(inputFileName);
117         setinp(inputFileName, nameLength, nameLength);
118         aread();
119         while(copen.nopen == -1) {
120                 copen.nopen=nopeno;
121             prnt1.iecho=1;
122                 xiniall();
123             aread();
124         }
125         Int_t ishini;
126         utpri("aamain",prnt1.ish,ishini,4,6);
127         if(appli.model != 1)
128                 IniModel(appli.model);
129         ebin.nrebin = 1;
130         ainit();
131         aseed(2);
132 }
133
134 void TEpos::GenerateEvent() {
135 //      cseed.seedj = gRandom->Rndm() * 1e10;
136         Int_t n = 1;
137         evgen(n);
138 }
139
140 Int_t TEpos::ImportParticles(TClonesArray *particles, Option_t *) {
141   //Fills provided ClonesArray with generated particles
142         particles->Clear();
143         if (!cevt.nevt) return 0;
144         Int_t numpart = cptl.nptl;
145         printf("%d particles generated\n", numpart);
146         for (Int_t iPart=0; iPart<numpart; iPart++) {
147                 Int_t tFather = cptl.iorptl[iPart] - 1;
148                 tFather = tFather < -1 ? -1 : tFather;
149                 if (tFather> -1) {
150                         TParticle *mother = (TParticle*) (particles->UncheckedAt(tFather));
151                         mother->SetLastDaughter(iPart);
152                         if (mother->GetFirstDaughter()==-1)
153                                 mother->SetFirstDaughter(iPart);
154                 }
155                 Int_t status = cptl.istptl[iPart] + 1;
156                 Int_t pdg = fIdConverter->ConvertIsajetToPdg(cptl.idptl[iPart]);
157                 new((*particles)[iPart]) TParticle(pdg, status,
158                                  tFather, -1, -1, -1,
159                                  cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
160                                  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);
161                 (*particles)[iPart]->SetUniqueID(iPart);
162         }
163
164         return numpart;
165 }
166
167 TObjArray*  TEpos::ImportParticles(Option_t *) {
168   //Creates new particle array
169         fParticles->Clear();
170         if (!cevt.nevt) return NULL;
171         Int_t numpart = cptl.nptl;
172         printf("%d particles generated\n", numpart);
173         for (Int_t iPart=0; iPart<numpart; iPart++) {
174                 Int_t tFather = cptl.iorptl[iPart] - 1;
175                 tFather = tFather < -1 ? -1 : tFather;
176                 if (tFather> -1) {
177                         TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
178                         mother->SetLastDaughter(iPart);
179                         if (mother->GetFirstDaughter()==-1)
180                                 mother->SetFirstDaughter(iPart);
181                 }
182                 Int_t status = cptl.istptl[iPart] + 1;
183                 Int_t pdg = fIdConverter->ConvertIsajetToPdg(cptl.idptl[iPart]);
184                 TParticle* p = new TParticle(pdg, status,
185                                  tFather, -1, -1, -1,
186                                  cptl.pptl[iPart][0], cptl.pptl[iPart][1],cptl.pptl[iPart][2],cptl.pptl[iPart][3],
187                                  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);
188                 p->SetUniqueID(iPart);
189                 fParticles->Add(p);
190         }
191
192         return fParticles;
193 }
194
195 void TEpos::AddNoDecay(Int_t nodecay) {
196         fNoDecays.push_back(nodecay);
197 }
198
199 void TEpos::AddExtraInputLine(const char *line) {
200         fExtraInputLines.push_back(line);
201 }
202
203 void TEpos::GenerateInputFile() {
204   // Generate input file in EPOS format
205         ofstream file(GetInputFileName(), ios_base::out | ios_base::trunc);
206         char epo[256];
207         char *epoEnv = getenv("EPO");
208         if (epoEnv) {
209                 strncpy(epo, epoEnv, 255);
210         } else {
211                 strncpy(epo, getenv("ALICE_ROOT"), 255);
212         }
213         strncat(epo, "/EPOS/epos167", 255);
214
215         file << "fname pathnx " << epo << "/" << endl;
216         file << "fname histo none" << endl;
217         file << "fname copy none" << endl;
218         file << "fname log none" << endl;
219         file << "fname check none" << endl;
220 //      file << "fname data /tmp/epos.out" << endl;
221         file << "fname initl " << epo << "/epos.initl" << endl;
222         file << "fname inidi " << epo << "/epos.inidi" << endl;
223         file << "fname inidr " << epo << "/epos.inidr" << endl;
224         file << "fname iniev " << epo << "/epos.iniev" << endl;
225         file << "fname inirj " << epo << "/epos.inirj" << endl;
226         file << "fname inics " << epo << "/epos.inics" << endl;
227         file << "fname inigrv " << epo << "/epos.inigrv" << endl;
228         file << "fqgsjet dat " << epo << "/qgsjet/qgsjet.dat" << endl;
229         file << "fqgsjet ncs " << epo << "/qgsjet/qgsjet.ncs" << endl;
230         file << "fqgsjetII dat " << epo << "/qgsjetII/qgsdat-II-03" << endl;
231         file << "fqgsjetII ncs " << epo << "/qgsjetII/sectnu-II-03" << endl;
232         file << "nodecay 120" << endl;
233         file << "nodecay -120" << endl;
234         file << "nodecay 130" << endl;
235         file << "nodecay -130" << endl;
236         file << "nodecay -20" << endl;
237         file << "nodecay 20" << endl;
238         file << "nodecay 14" << endl;
239         file << "nodecay -14" << endl;
240         file << "set ndecay 1111110" << endl;
241         file << "echo on" << endl;
242
243 //      .optns file
244         int precision = file.precision();
245         file.precision(15);
246         file << "set seedj " << (gRandom->Rndm() * 1e14) << endl;
247         file.precision(precision);
248         file << "application hadron" << endl;
249         file << "set laproj " << fLaproj << endl;
250         file << "set maproj " << fMaproj << endl;
251         file << "set latarg " << fLatarg << endl;
252         file << "set matarg " << fMatarg << endl;
253         file << "set bminim " << fBminim << endl;
254         file << "set bmaxim " << fBmaxim << endl;
255         file << "set phimin " << fPhimin << endl;
256         file << "set phimax " << fPhimax << endl;
257         file << "set ecms " << fEcms << endl;
258
259         for(unsigned int i = 0; i < fNoDecays.size(); ++i) {
260                 file << "nodecay " << fNoDecays[i] << endl;
261         }
262
263         file << "switch splitting " << (fSplitting ? "on" : "off") << endl;
264
265         file << "frame nucleon-nucleon" << endl;
266
267         for(unsigned int i = 0; i < fExtraInputLines.size(); ++i) {
268                 file << fExtraInputLines[i] << endl;
269         }
270
271 //    file << "output epos" << endl;
272 //    file << "record event nevt nptl b endrecord" << endl;
273 //    file << "record particle i id fa mo c1 c2 st endrecord" << endl;
274
275         file << "input " << epo << "/epos.param" << endl;
276         file << "runprogram" << endl;
277         file.close();
278 }
279
280 Float_t TEpos::GetBimevt() const { return cevt.bimevt; }
281 Float_t TEpos::GetPhievt() const { return cevt.phievt; }
282 Int_t TEpos::GetKolevt() const { return cevt.kolevt; }
283 Int_t TEpos::GetKoievt() const { return cevt.koievt; }
284 Float_t TEpos::GetPmxevt() const { return cevt.pmxevt; }
285 Float_t TEpos::GetEgyevt() const { return cevt.egyevt; }
286 Int_t TEpos::GetNpjevt() const { return cevt.npjevt; }
287 Int_t TEpos::GetNtgevt() const { return cevt.ntgevt; }
288 Int_t TEpos::GetNpnevt() const { return cevt.npnevt; }
289 Int_t TEpos::GetNppevt() const { return cevt.nppevt; }
290 Int_t TEpos::GetNtnevt() const { return cevt.ntnevt; }
291 Int_t TEpos::GetNtpevt() const { return cevt.ntpevt; }
292 Int_t TEpos::GetJpnevt() const { return cevt.jpnevt; }
293 Int_t TEpos::GetJppevt() const { return cevt.jppevt; }
294 Int_t TEpos::GetJtnevt() const { return cevt.jtnevt; }
295 Int_t TEpos::GetJtpevt() const { return cevt.jtpevt; }
296 Float_t TEpos::GetXbjevt() const { return cevt.xbjevt; }
297 Float_t TEpos::GetQsqevt() const { return cevt.qsqevt; }
298 Int_t TEpos::GetNglevt() const { return cevt.nglevt; }
299 Float_t TEpos::GetZppevt() const { return cevt.zppevt; }
300 Float_t TEpos::GetZptevt() const { return cevt.zptevt; }
301