TTherminator/AliGenTherminator.cxx

   1 // ALICE event generator based on the THERMINATOR model
   2 // It reads the test output of the model and puts it onto
   3 // the stack
   4 // It has an option to use the Lhyquid3D input freeze-out
   5 // hypersurface
   6 // Author: Adam.Kisiel@cern.ch
   7
   8 #include <iostream>
   9 #include <fstream>
  10 #include <sstream>
  11 #include <TClonesArray.h>
  12 #include <TMCProcess.h>
  13 #include <TDatabasePDG.h>
  14 #include <TParticle.h>
  15
  16 #include "AliConst.h"
  17 #include "AliDecayer.h"
  18 #include "AliGenEventHeader.h"
  19 #include "AliGenHijingEventHeader.h"
  20 #include "AliGenTherminator.h"
  21 #include "AliLog.h"
  22 #include "AliRun.h"
  23
  24 ClassImp(AliGenTherminator)
  25
  26 using namespace std;
  27
  28 AliGenTherminator::AliGenTherminator():
  29   AliGenMC(),
  30   fNt(0),
  31   fEventNumber(0),
  32   fFileName(""),
  33   fFreezeOutModel(""),
  34   fTemperature(0.1656),
  35   fMiuI(-0.0009),
  36   fMiuS(0.0),
  37   fMiuB(0.0008),
  38   fAlfaRange(6.0),
  39   fRapRange(4.0),
  40   fRhoMax(8.0),
  41   fTau(8.0),
  42   fBWA(0.0),
  43   fBWVt(1.41),
  44   fBWDelay(0.0)
  45 {
  46   // Default constructor
  47   fParticles = new TClonesArray("TParticle",20000);
  48 }
  49 AliGenTherminator::AliGenTherminator(Int_t npart):
  50   AliGenMC(npart),
  51   fNt(0),
  52   fEventNumber(0),
  53   fFileName(""),
  54   fFreezeOutModel(""),
  55   fTemperature(0.1656),
  56   fMiuI(-0.0009),
  57   fMiuS(0.0),
  58   fMiuB(0.0008),
  59   fAlfaRange(6.0),
  60   fRapRange(4.0),
  61   fRhoMax(8.0),
  62   fTau(8.0),
  63   fBWA(0.0),
  64   fBWVt(1.41),
  65   fBWDelay(0.0)
  66 {
  67   // Constructor specifying the size of the particle table
  68   fParticles = new TClonesArray("TParticle",20000);
  69   fNprimaries = 0;
  70 }
  71
  72 AliGenTherminator::~AliGenTherminator()
  73 {
  74   //  AliGenMC::~AliGenMC();
  75   //  if (fTherminator) delete fTherminator;
  76 }
  77
  78 void AliGenTherminator::Generate()
  79 {
  80   // Run single event generation with the Therminator model
  81   AliWarning("Generating event from AliGenTherminator");
  82
  83   Float_t polar[3]    =   {0,0,0};
  84   Float_t origin[3]   =   {0,0,0};
  85   Float_t origin0[3]  =   {0,0,0};
  86   Float_t p[3];
  87   Float_t mass, energy;
  88
  89   Int_t nt  = 0;
  90   Int_t j, kf, ks, imo;
  91   kf = 0;
  92
  93   Vertex();
  94   for (j=0; j < 3; j++) origin0[j] = fVertex[j];
  95
  96   // Generate one event
  97
  98   ((TTherminator *) fMCEvGen)->GenerateEvent();
  99   AliWarning("Generated");
 100   ((TTherminator *) fMCEvGen)->ImportParticles(fParticles);
 101
 102   Int_t np = fParticles->GetEntriesFast();
 103   AliWarning(Form("Imported %d particles", np));
 104
 105   TParticle *iparticle;
 106   Double_t evrot = gRandom->Rndm()*TMath::Pi();
 107
 108   for (int i = 0; i < np; i++) {
 109     iparticle = (TParticle *) fParticles->At(i);
 110     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
 111     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
 112
 113     kf   = iparticle->GetPdgCode();
 114     ks   = iparticle->GetStatusCode();
 115     Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
 116     Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
 117     p[0] = arho*TMath::Cos(aphi + evrot);
 118     p[1] = arho*TMath::Sin(aphi + evrot);
 119 //     p[0] = iparticle->Px();
 120 //     p[1] = iparticle->Py();
 121     p[2] = iparticle->Pz();
 122     mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
 123     energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 124     origin[0] = origin0[0]+iparticle->Vx();
 125     origin[1] = origin0[1]+iparticle->Vy();
 126     origin[2] = origin0[2]+iparticle->Vz();
 127
 128     imo = -1;
 129     TParticle* mother = 0;
 130     if (hasMother) {
 131       imo = iparticle->GetFirstMother();
 132       mother = (TParticle *) fParticles->At(imo);
 133     } // if has mother
 134     Bool_t tFlag = (!hasDaughter);
 135
 136     printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
 137     PushTrack(tFlag,imo,kf,
 138               p[0],p[1],p[2],energy,
 139               origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
 140               polar[0],polar[1],polar[2],
 141               hasMother ? kPDecay:kPNoProcess,nt);
 142
 143     fNprimaries++;
 144     KeepTrack(nt);
 145   }
 146
 147   SetHighWaterMark(fNprimaries);
 148
 149   TArrayF eventVertex;
 150   eventVertex.Set(3);
 151   eventVertex[0] = origin0[0];
 152   eventVertex[1] = origin0[1];
 153   eventVertex[2] = origin0[2];
 154
 155 // Builds the event header, to be called after each event
 156   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
 157
 158   // Header
 159   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
 160   // Event Vertex
 161 //   header->SetPrimaryVertex(eventVertex);
 162 //   header->SetNProduced(fNprimaries);
 163
 164   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
 165   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
 166   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
 167   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
 168   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
 169   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
 170   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
 171   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
 172   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
 173
 174
 175 // 4-momentum vectors of the triggered jets.
 176 //
 177 // Before final state gluon radiation.
 178 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
 179 //                                            fHijing->GetHINT1(22),
 180 //                                            fHijing->GetHINT1(23),
 181 //                                            fHijing->GetHINT1(24));
 182
 183 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
 184 //                                            fHijing->GetHINT1(32),
 185 //                                            fHijing->GetHINT1(33),
 186 //                                            fHijing->GetHINT1(34));
 187 // // After final state gluon radiation.
 188 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
 189 //                                            fHijing->GetHINT1(27),
 190 //                                            fHijing->GetHINT1(28),
 191 //                                            fHijing->GetHINT1(29));
 192
 193 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
 194 //                                            fHijing->GetHINT1(37),
 195 //                                            fHijing->GetHINT1(38),
 196 //                                            fHijing->GetHINT1(39));
 197 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
 198 // Bookkeeping for kinematic bias
 199 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
 200 // Event Vertex
 201   header->SetPrimaryVertex(fVertex);
 202   AddHeader(header);
 203   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 204
 205   gAlice->SetGenEventHeader(header);
 206 }
 207
 208 void AliGenTherminator::Init()
 209 {
 210   // Initialize global variables and
 211   // particle and decay tables
 212   if (fFileName.Length() == 0)
 213     fFileName = "event.out";
 214   ReadShareParticleTable();
 215
 216   SetMC(new TTherminator());
 217
 218   AliGenMC::Init();
 219   ((TTherminator *) fMCEvGen)->Initialize();
 220 }
 221 void AliGenTherminator::SetFileName(const char *infilename)
 222 {
 223   // Set parameter filename
 224   fFileName = infilename;
 225 }
 226 void AliGenTherminator::SetEventNumberInFile(int evnum)
 227 {
 228   // Set number of events to generate - default: 1
 229   fEventNumber = evnum;
 230 }
 231
 232 void AliGenTherminator::ReadShareParticleTable()
 233 {
 234   // Read in particle table from share
 235   // and add missing particle type to TDatabasePDG
 236
 237   char str[50];
 238   char str1[200];
 239
 240   TDatabasePDG *tInstance = TDatabasePDG::Instance();
 241   TParticlePDG *tParticleType;
 242
 243   AliWarning(Form("Reading particle types from particles.data"));
 244
 245   TString aroot = gSystem->Getenv("ALICE_ROOT");
 246   ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
 247   //  ifstream in("particles.data");
 248
 249   int charge;
 250
 251   int number=0;
 252   if ((in) && (in.is_open()))
 253     {
 254       //START OF HEAD-LINE
 255       in.ignore(200,'\n');
 256       in.ignore(200,'\n');
 257       in.ignore(200,'\n');
 258       //END OF HEAD-LINE
 259
 260       while (in>>str)
 261         {
 262           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 263             {
 264               in.getline(str1,200);
 265               continue;
 266             }
 267           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
 268
 269           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
 270           number++;
 271           tParticleType = tInstance->GetParticle((int) mc);
 272           if (!tParticleType) {
 273             charge = 0;
 274             if (strstr(str, "plu")) charge = 1;
 275             if (strstr(str, "min")) charge = -1;
 276             if (strstr(str, "plb")) charge = -1;
 277             if (strstr(str, "mnb")) charge = 1;
 278             if (strstr(str, "plp")) charge = 2;
 279             if (strstr(str, "ppb")) charge = -2;
 280             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
 281             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
 282             //      AliWarning(Form("Quantum numbers q s c aq as ac tI3 %lf %lf %lf %lf %lf %lf %lf", q, s, c, aq, as, ac, tI3));
 283
 284           }
 285         }
 286       in.close();
 287     }
 288   CreateTherminatorInputFile();
 289 }
 290
 291 void AliGenTherminator::CreateTherminatorInputFile()
 292 {
 293   // Create Therminator input file
 294   const char *aroot = gSystem->Getenv("ALICE_ROOT");
 295   ofstream *ostr = new ofstream("therminator.in");
 296   (*ostr) << "NumberOfEvents = 1" << endl;
 297   (*ostr) << "Randomize = 1" << endl;
 298   (*ostr) << "TableType = SHARE" << endl;
 299   (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
 300   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
 301   (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
 302   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
 303   (*ostr) << "BWVt = " << fBWVt << endl;
 304   (*ostr) << "Tau = " << fTau << endl;
 305   (*ostr) << "RhoMax = " << fRhoMax << endl;
 306   (*ostr) << "Temperature = " << fTemperature << endl;
 307   (*ostr) << "MiuI = " << fMiuI << endl;
 308   (*ostr) << "MiuS = " << fMiuS << endl;
 309   (*ostr) << "MiuB = " << fMiuB << endl;
 310   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
 311   (*ostr) << "RapidityRange = " << fRapRange << endl;
 312   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
 313 }
 314
 315 void AliGenTherminator::SetModel(const char *model)
 316 {
 317   // Set the freeze-out model to use
 318   fFreezeOutModel = model;
 319 }
 320
 321 void AliGenTherminator::SetLhyquidSet(const char *set)
 322 {
 323   // Select one of pregenerated Lhyquid hypersurfaces
 324   const char *aroot = gSystem->Getenv("ALICE_ROOT");
 325   if (strstr(set, "LHC500C0005")) {
 326     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
 327     AliWarning(Form("  Pb-Pb collisions, centrality 0-5%"));
 328     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
 329     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
 330     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
 331     fFOHSlocation = aroot;
 332     fFOHSlocation += "/TTherminator/data/LHC500C0005";
 333   }
 334   else if (strstr(set, "LHC500C2030")) {
 335     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
 336     AliWarning(Form("  Pb-Pb collisions, centrality 20-30%"));
 337     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
 338     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
 339     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
 340     fFOHSlocation = aroot;
 341     fFOHSlocation += "/TTherminator/data/LHC500C2030";
 342   }
 343   else {
 344     AliWarning(Form("Did not find Lhyquid set %s", set));
 345     AliWarning(Form("Reverting to default: current directory"));
 346     fFOHSlocation = "";
 347   }
 348 }
 349
 350 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
 351 {
 352   // Select Your own Lhyquid hypersurface
 353   fFOHSlocation = inputdir;
 354 }