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 // Author: Adam.Kisiel@cern.ch
   5
   6 #include <iostream>
   7 #include <fstream>
   8 #include <sstream>
   9 #include <TClonesArray.h>
  10 #include <TMCProcess.h>
  11 #include <TDatabasePDG.h>
  12 #include <TParticle.h>
  13
  14 #include "AliConst.h"
  15 #include "AliDecayer.h"
  16 #include "AliGenEventHeader.h"
  17 #include "AliGenHijingEventHeader.h"
  18 #include "AliGenTherminator.h"
  19 #include "AliLog.h"
  20 #include "AliRun.h"
  21
  22 ClassImp(AliGenTherminator)
  23
  24 using namespace std;
  25
  26 AliGenTherminator::AliGenTherminator():
  27   AliGenMC(),
  28   fNt(0),
  29   fEventNumber(0),
  30   fFileName(""),
  31   fFreezeOutModel(""),
  32   fTemperature(0.1656),
  33   fMiuI(-0.0009),
  34   fMiuS(0.0),
  35   fMiuB(0.0008),
  36   fAlfaRange(6.0),
  37   fRapRange(4.0),
  38   fRhoMax(8.0),
  39   fTau(8.0),
  40   fBWA(0.0),
  41   fBWVt(1.41),
  42   fBWDelay(0.0)
  43 {
  44   // Default constructor
  45   fParticles = new TClonesArray("TParticle",20000);
  46 }
  47 AliGenTherminator::AliGenTherminator(Int_t npart):
  48   AliGenMC(npart),
  49   fNt(0),
  50   fEventNumber(0),
  51   fFileName(""),
  52   fFreezeOutModel(""),
  53   fTemperature(0.1656),
  54   fMiuI(-0.0009),
  55   fMiuS(0.0),
  56   fMiuB(0.0008),
  57   fAlfaRange(6.0),
  58   fRapRange(4.0),
  59   fRhoMax(8.0),
  60   fTau(8.0),
  61   fBWA(0.0),
  62   fBWVt(1.41),
  63   fBWDelay(0.0)
  64 {
  65   // Constructor specifying the size of the particle table
  66   fParticles = new TClonesArray("TParticle",20000);
  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   TParticle *iparticle;
 104   Double_t evrot = gRandom->Rndm()*TMath::Pi();
 105
 106   for (int i = 0; i < np; i++) {
 107     iparticle = (TParticle *) fParticles->At(i);
 108     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
 109     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
 110
 111     kf   = iparticle->GetPdgCode();
 112     ks   = iparticle->GetStatusCode();
 113     Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
 114     Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
 115     p[0] = arho*TMath::Cos(aphi + evrot);
 116     p[1] = arho*TMath::Sin(aphi + evrot);
 117 //     p[0] = iparticle->Px();
 118 //     p[1] = iparticle->Py();
 119     p[2] = iparticle->Pz();
 120     mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
 121     energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 122     origin[0] = origin0[0]+iparticle->Vx();
 123     origin[1] = origin0[1]+iparticle->Vy();
 124     origin[2] = origin0[2]+iparticle->Vz();
 125
 126     imo = -1;
 127     TParticle* mother = 0;
 128     if (hasMother) {
 129       imo = iparticle->GetFirstMother();
 130       mother = (TParticle *) fParticles->At(imo);
 131     } // if has mother
 132     Bool_t tFlag = (!hasDaughter);
 133
 134     printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
 135     PushTrack(tFlag,imo,kf,
 136               p[0],p[1],p[2],energy,
 137               origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
 138               polar[0],polar[1],polar[2],
 139               hasMother ? kPDecay:kPNoProcess,nt);
 140
 141     fNprimaries++;
 142     KeepTrack(nt);
 143   }
 144
 145   SetHighWaterMark(fNprimaries);
 146
 147   TArrayF eventVertex;
 148   eventVertex.Set(3);
 149   eventVertex[0] = origin0[0];
 150   eventVertex[1] = origin0[1];
 151   eventVertex[2] = origin0[2];
 152
 153 // Builds the event header, to be called after each event
 154   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
 155
 156   // Header
 157   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
 158   // Event Vertex
 159 //   header->SetPrimaryVertex(eventVertex);
 160 //   header->SetNProduced(fNprimaries);
 161
 162   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
 163   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
 164   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
 165   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
 166   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
 167   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
 168   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
 169   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
 170   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
 171
 172
 173 // 4-momentum vectors of the triggered jets.
 174 //
 175 // Before final state gluon radiation.
 176 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
 177 //                                            fHijing->GetHINT1(22),
 178 //                                            fHijing->GetHINT1(23),
 179 //                                            fHijing->GetHINT1(24));
 180
 181 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
 182 //                                            fHijing->GetHINT1(32),
 183 //                                            fHijing->GetHINT1(33),
 184 //                                            fHijing->GetHINT1(34));
 185 // // After final state gluon radiation.
 186 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
 187 //                                            fHijing->GetHINT1(27),
 188 //                                            fHijing->GetHINT1(28),
 189 //                                            fHijing->GetHINT1(29));
 190
 191 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
 192 //                                            fHijing->GetHINT1(37),
 193 //                                            fHijing->GetHINT1(38),
 194 //                                            fHijing->GetHINT1(39));
 195 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
 196 // Bookkeeping for kinematic bias
 197 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
 198 // Event Vertex
 199   header->SetPrimaryVertex(fVertex);
 200   AddHeader(header);
 201   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 202
 203   gAlice->SetGenEventHeader(header);
 204 }
 205
 206 void AliGenTherminator::Init()
 207 {
 208   // Initialize global variables and
 209   // particle and decay tables
 210   if (fFileName.Length() == 0)
 211     fFileName = "event.out";
 212   ReadShareParticleTable();
 213
 214   SetMC(new TTherminator());
 215
 216   AliGenMC::Init();
 217   ((TTherminator *) fMCEvGen)->Initialize();
 218 }
 219 void AliGenTherminator::SetFileName(const char *infilename)
 220 {
 221   // Set parameter filename
 222   fFileName = infilename;
 223 }
 224 void AliGenTherminator::SetEventNumberInFile(int evnum)
 225 {
 226   // Set number of events to generate - default: 1
 227   fEventNumber = evnum;
 228 }
 229
 230 void AliGenTherminator::ReadShareParticleTable()
 231 {
 232   // Read in particle table from share
 233   // and add missing particle type to TDatabasePDG
 234
 235   char str[50];
 236   char str1[200];
 237
 238   TDatabasePDG *tInstance = TDatabasePDG::Instance();
 239   TParticlePDG *tParticleType;
 240
 241   AliWarning(Form("Reading particle types from particles.data"));
 242   ifstream in("particles.data");
 243
 244   int charge;
 245
 246   int number=0;
 247   if ((in) && (in.is_open()))
 248     {
 249       //START OF HEAD-LINE
 250       in.ignore(200,'\n');
 251       in.ignore(200,'\n');
 252       in.ignore(200,'\n');
 253       //END OF HEAD-LINE
 254
 255       while (in>>str)
 256         {
 257           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 258             {
 259               in.getline(str1,200);
 260               continue;
 261             }
 262           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
 263
 264           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
 265           number++;
 266           tParticleType = tInstance->GetParticle((int) mc);
 267           if (!tParticleType) {
 268             charge = 0;
 269             if (strstr(str, "plu")) charge = 1;
 270             if (strstr(str, "min")) charge = -1;
 271             if (strstr(str, "plb")) charge = -1;
 272             if (strstr(str, "mnb")) charge = 1;
 273             if (strstr(str, "plp")) charge = 2;
 274             if (strstr(str, "ppb")) charge = -2;
 275             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
 276             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
 277             //      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));
 278
 279           }
 280         }
 281       in.close();
 282     }
 283   CreateTherminatorInputFile();
 284 }
 285
 286 void AliGenTherminator::CreateTherminatorInputFile()
 287 {
 288   // Create Therminator input file
 289   ofstream *ostr = new ofstream("therminator.in");
 290   (*ostr) << "NumberOfEvents = 1" << endl;
 291   (*ostr) << "Randomize = 1" << endl;
 292   (*ostr) << "TableType = SHARE" << endl;
 293   (*ostr) << "InputDirSHARE = ." << endl;
 294   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
 295   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
 296   (*ostr) << "BWVt = " << fBWVt << endl;
 297   (*ostr) << "Tau = " << fTau << endl;
 298   (*ostr) << "RhoMax = " << fRhoMax << endl;
 299   (*ostr) << "Temperature = " << fTemperature << endl;
 300   (*ostr) << "MiuI = " << fMiuI << endl;
 301   (*ostr) << "MiuS = " << fMiuS << endl;
 302   (*ostr) << "MiuB = " << fMiuB << endl;
 303   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
 304   (*ostr) << "RapidityRange = " << fRapRange << endl;
 305   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
 306 }
 307
 308 void AliGenTherminator::SetModel(const char *model)
 309 {
 310   // Set the freeze-out model to use
 311   fFreezeOutModel = model;
 312 }