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