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 }
  48 AliGenTherminator::AliGenTherminator(Int_t npart):
  49   AliGenMC(npart),
  50   fNt(0),
  51   fEventNumber(0),
  52   fFileName(""),
  53   fFreezeOutModel(""),
  54   fTemperature(0.1656),
  55   fMiuI(-0.0009),
  56   fMiuS(0.0),
  57   fMiuB(0.0008),
  58   fAlfaRange(6.0),
  59   fRapRange(4.0),
  60   fRhoMax(8.0),
  61   fTau(8.0),
  62   fBWA(0.0),
  63   fBWVt(1.41),
  64   fBWDelay(0.0)
  65 {
  66   // Constructor specifying the size of the particle table
  67   fNprimaries = 0;
  68 }
  69
  70 AliGenTherminator::~AliGenTherminator()
  71 {
  72   //  AliGenMC::~AliGenMC();
  73   //  if (fTherminator) delete fTherminator;
  74 }
  75
  76 void AliGenTherminator::Generate()
  77 {
  78   // Run single event generation with the Therminator model
  79   AliWarning("Generating event from AliGenTherminator");
  80
  81   Float_t polar[3]    =   {0,0,0};
  82   Float_t origin[3]   =   {0,0,0};
  83   Float_t origin0[3]  =   {0,0,0};
  84   Float_t p[3];
  85   Float_t mass, energy;
  86
  87   Int_t nt  = 0;
  88   Int_t j, kf, ks, imo;
  89   kf = 0;
  90
  91   Vertex();
  92   for (j=0; j < 3; j++) origin0[j] = fVertex[j];
  93
  94   // Generate one event
  95
  96   ((TTherminator *) fMCEvGen)->GenerateEvent();
  97   AliWarning("Generated");
  98   ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
  99
 100   Int_t np = fParticles.GetEntriesFast();
 101   AliWarning(Form("Imported %d particles", np));
 102
 103   Int_t *idsOnStack;
 104   idsOnStack = new Int_t[np];
 105
 106   TParticle *iparticle;
 107   Double_t evrot = gRandom->Rndm()*TMath::Pi();
 108
 109   for (int i = 0; i < np; i++) {
 110     iparticle = (TParticle *) fParticles.At(i);
 111     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
 112     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
 113
 114     if (hasDaughter) {
 115       // This particle has decayed
 116       // It will not be tracked
 117       // Add it only once with coorduinates not
 118       // smeared with primary vertex position
 119
 120       kf   = iparticle->GetPdgCode();
 121       ks   = iparticle->GetStatusCode();
 122       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
 123       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
 124       p[0] = arho*TMath::Cos(aphi + evrot);
 125       p[1] = arho*TMath::Sin(aphi + evrot);
 126       //     p[0] = iparticle->Px();
 127       //     p[1] = iparticle->Py();
 128       p[2] = iparticle->Pz();
 129       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
 130       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 131
 132       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
 133       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
 134       origin[0] = vrho*TMath::Cos(vphi + evrot);
 135       origin[1] = vrho*TMath::Sin(vphi + evrot);
 136       origin[2] = iparticle->Vz();
 137
 138       imo = -1;
 139       TParticle* mother = 0;
 140       if (hasMother) {
 141         imo = iparticle->GetFirstMother();
 142         mother = (TParticle *) fParticles.At(imo);
 143       } // if has mother
 144       Bool_t tFlag = (!hasDaughter);
 145
 146       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
 147       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
 148                 p[0],p[1],p[2],energy,
 149                 origin[0],origin[1],origin[2],iparticle->T(),
 150                 polar[0],polar[1],polar[2],
 151                 hasMother ? kPDecay:kPNoProcess,nt);
 152       idsOnStack[i] = nt;
 153       fNprimaries++;
 154       KeepTrack(nt);
 155     }
 156     else {
 157       // This is a final state particle
 158       // It will be tracked
 159       // Add it TWICE to the stack !!!
 160       // First time with event-wide coordicates (for femtoscopy) -
 161       //   this one will not be tracked
 162       // Second time with event-wide ccordiantes and vertex smearing
 163       //   this one will be tracked
 164
 165       kf   = iparticle->GetPdgCode();
 166       ks   = iparticle->GetStatusCode();
 167       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
 168       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
 169       p[0] = arho*TMath::Cos(aphi + evrot);
 170       p[1] = arho*TMath::Sin(aphi + evrot);
 171       //     p[0] = iparticle->Px();
 172       //     p[1] = iparticle->Py();
 173       p[2] = iparticle->Pz();
 174       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
 175       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 176
 177       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
 178       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
 179       origin[0] = vrho*TMath::Cos(vphi + evrot);
 180       origin[1] = vrho*TMath::Sin(vphi + evrot);
 181       origin[2] = iparticle->Vz();
 182
 183       imo = -1;
 184       TParticle* mother = 0;
 185       if (hasMother) {
 186         imo = iparticle->GetFirstMother();
 187         mother = (TParticle *) fParticles.At(imo);
 188       } // if has mother
 189       Bool_t tFlag = (hasDaughter);
 190
 191       printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
 192       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
 193       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
 194                 p[0],p[1],p[2],energy,
 195                 origin[0],origin[1],origin[2],iparticle->T(),
 196                 polar[0],polar[1],polar[2],
 197                 hasMother ? kPDecay:kPNoProcess,nt);
 198       idsOnStack[i] = nt;
 199       fNprimaries++;
 200       KeepTrack(nt);
 201
 202       origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
 203       origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
 204       origin[2] = origin0[2]+iparticle->Vz();
 205
 206       imo = nt;
 207       //      mother = (TParticle *) fParticles.At(nt);
 208       tFlag = (!hasDaughter);
 209
 210       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
 211       PushTrack(tFlag,imo,kf,
 212                 p[0],p[1],p[2],energy,
 213                 origin[0],origin[1],origin[2],iparticle->T(),
 214                 polar[0],polar[1],polar[2],
 215                 hasMother ? kPDecay:kPNoProcess,nt);
 216       fNprimaries++;
 217       KeepTrack(nt);
 218     }
 219   }
 220
 221
 222
 223   SetHighWaterMark(fNprimaries);
 224
 225   TArrayF eventVertex;
 226   eventVertex.Set(3);
 227   eventVertex[0] = origin0[0];
 228   eventVertex[1] = origin0[1];
 229   eventVertex[2] = origin0[2];
 230
 231 // Builds the event header, to be called after each event
 232   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
 233
 234   // Header
 235   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
 236   // Event Vertex
 237 //   header->SetPrimaryVertex(eventVertex);
 238 //   header->SetNProduced(fNprimaries);
 239
 240   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
 241   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
 242   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
 243   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
 244   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
 245   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
 246   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
 247   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
 248   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
 249
 250
 251 // 4-momentum vectors of the triggered jets.
 252 //
 253 // Before final state gluon radiation.
 254 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
 255 //                                            fHijing->GetHINT1(22),
 256 //                                            fHijing->GetHINT1(23),
 257 //                                            fHijing->GetHINT1(24));
 258
 259 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
 260 //                                            fHijing->GetHINT1(32),
 261 //                                            fHijing->GetHINT1(33),
 262 //                                            fHijing->GetHINT1(34));
 263 // // After final state gluon radiation.
 264 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
 265 //                                            fHijing->GetHINT1(27),
 266 //                                            fHijing->GetHINT1(28),
 267 //                                            fHijing->GetHINT1(29));
 268
 269 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
 270 //                                            fHijing->GetHINT1(37),
 271 //                                            fHijing->GetHINT1(38),
 272 //                                            fHijing->GetHINT1(39));
 273 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
 274 // Bookkeeping for kinematic bias
 275 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
 276 // Event Vertex
 277   header->SetPrimaryVertex(fVertex);
 278   AddHeader(header);
 279   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 280
 281   delete idsOnStack;
 282
 283   //  gAlice->SetGenEventHeader(header);
 284 }
 285
 286 void AliGenTherminator::Init()
 287 {
 288   // Initialize global variables and
 289   // particle and decay tables
 290   if (fFileName.Length() == 0)
 291     fFileName = "event.out";
 292   ReadShareParticleTable();
 293
 294   SetMC(new TTherminator());
 295
 296   AliGenMC::Init();
 297   ((TTherminator *) fMCEvGen)->Initialize();
 298 }
 299 void AliGenTherminator::SetFileName(const char *infilename)
 300 {
 301   // Set parameter filename
 302   fFileName = infilename;
 303 }
 304 void AliGenTherminator::SetEventNumberInFile(int evnum)
 305 {
 306   // Set number of events to generate - default: 1
 307   fEventNumber = evnum;
 308 }
 309
 310 void AliGenTherminator::ReadShareParticleTable()
 311 {
 312   // Read in particle table from share
 313   // and add missing particle type to TDatabasePDG
 314
 315   char str[50];
 316   char str1[200];
 317
 318   TDatabasePDG *tInstance = TDatabasePDG::Instance();
 319   TParticlePDG *tParticleType;
 320
 321   AliWarning(Form("Reading particle types from particles.data"));
 322
 323   TString aroot = gSystem->Getenv("ALICE_ROOT");
 324   ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
 325   //  ifstream in("particles.data");
 326
 327   int charge;
 328
 329   int number=0;
 330   if ((in) && (in.is_open()))
 331     {
 332       //START OF HEAD-LINE
 333       in.ignore(200,'\n');
 334       in.ignore(200,'\n');
 335       in.ignore(200,'\n');
 336       //END OF HEAD-LINE
 337
 338       while (in>>str)
 339         {
 340           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 341             {
 342               in.getline(str1,200);
 343               continue;
 344             }
 345           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
 346
 347           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
 348           number++;
 349           tParticleType = tInstance->GetParticle((int) mc);
 350           if (!tParticleType) {
 351             charge = 0;
 352             if (strstr(str, "plu")) charge = 1;
 353             if (strstr(str, "min")) charge = -1;
 354             if (strstr(str, "plb")) charge = -1;
 355             if (strstr(str, "mnb")) charge = 1;
 356             if (strstr(str, "plp")) charge = 2;
 357             if (strstr(str, "ppb")) charge = -2;
 358             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
 359             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
 360             //      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));
 361
 362           }
 363         }
 364       in.close();
 365     }
 366   CreateTherminatorInputFile();
 367 }
 368
 369 void AliGenTherminator::CreateTherminatorInputFile()
 370 {
 371   // Create Therminator input file
 372   const char *aroot = gSystem->Getenv("ALICE_ROOT");
 373   ofstream *ostr = new ofstream("therminator.in");
 374   (*ostr) << "NumberOfEvents = 1" << endl;
 375   (*ostr) << "Randomize = 1" << endl;
 376   (*ostr) << "TableType = SHARE" << endl;
 377   (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
 378   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
 379   (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
 380   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
 381   (*ostr) << "BWVt = " << fBWVt << endl;
 382   (*ostr) << "Tau = " << fTau << endl;
 383   (*ostr) << "RhoMax = " << fRhoMax << endl;
 384   (*ostr) << "Temperature = " << fTemperature << endl;
 385   (*ostr) << "MiuI = " << fMiuI << endl;
 386   (*ostr) << "MiuS = " << fMiuS << endl;
 387   (*ostr) << "MiuB = " << fMiuB << endl;
 388   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
 389   (*ostr) << "RapidityRange = " << fRapRange << endl;
 390   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
 391 }
 392
 393 void AliGenTherminator::SetModel(const char *model)
 394 {
 395   // Set the freeze-out model to use
 396   fFreezeOutModel = model;
 397 }
 398
 399 void AliGenTherminator::SetLhyquidSet(const char *set)
 400 {
 401   // Select one of pregenerated Lhyquid hypersurfaces
 402   const char *aroot = gSystem->Getenv("ALICE_ROOT");
 403   if (strstr(set, "LHC500C0005")) {
 404     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
 405     AliWarning(Form("  Pb-Pb collisions, centrality 0-5%"));
 406     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
 407     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
 408     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
 409     fFOHSlocation = aroot;
 410     fFOHSlocation += "/TTherminator/data/LHC500C0005";
 411   }
 412   else if (strstr(set, "LHC500C2030")) {
 413     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
 414     AliWarning(Form("  Pb-Pb collisions, centrality 20-30%"));
 415     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
 416     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
 417     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
 418     fFOHSlocation = aroot;
 419     fFOHSlocation += "/TTherminator/data/LHC500C2030";
 420   }
 421   else {
 422     AliWarning(Form("Did not find Lhyquid set %s", set));
 423     AliWarning(Form("Reverting to default: current directory"));
 424     fFOHSlocation = "";
 425   }
 426 }
 427
 428 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
 429 {
 430   // Select Your own Lhyquid hypersurface
 431   fFOHSlocation = inputdir;
 432 }