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