TTherminator/TTherminator.cxx

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //                                                                           //
   3 // TTherminator: global interface class to the Therminator model.            //
   4 // Initialized global variables and runs the event generation                //
   5 //                                                                           //
   6 ///////////////////////////////////////////////////////////////////////////////
   7 /******************************************************************************
   8  *                      T H E R M I N A T O R                                 *
   9  *                   THERMal heavy-IoN generATOR                              *
  10  *                           version 1.0                                      *
  11  *                                                                            *
  12  * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
  13  *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
  14  * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
  15  *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
  16  * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
  17  *                 Wojciech Florkowski                                        *
  18  *                                                                            *
  19  * For the detailed description of the program and furhter references         *
  20  * to the description of the model plesase refer to: nucl-th/0504047,         *
  21  * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
  22  *                                                                            *
  23  * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
  24  *                                                                            *
  25  * This code can be freely used and redistributed. However if you decide to   *
  26  * make modifications to the code, please contact the authors, especially     *
  27  * if you plan to publish the results obtained with such modified code.       *
  28  * Any publication of results obtained using this code must include the       *
  29  * reference to nucl-th/0504047 and the published version of it, when         *
  30  * available.                                                                 *
  31  *                                                                            *
  32  *****************************************************************************/
  33 #include <TTherminator.h>
  34 #include <fstream>
  35 #include <TDatabasePDG.h>
  36 #include <TParticle.h>
  37 #include <TClonesArray.h>
  38 #include <TClass.h>
  39
  40 ReadPar *sRPInstance;
  41 STR      sRPFileName;
  42 int      sNumEvents;
  43 int      sRunType;
  44 int      sTables;
  45 int      sModel;
  46 int      sIntegrateSample;
  47
  48 ClassImp(TTherminator)
  49
  50 TTherminator::TTherminator():
  51   fCalka(0),
  52   fEvent(0),
  53   fPartDB(0)
  54 {
  55   // Default constructor
  56   fPartDB = new ParticleDB();
  57 }
  58 TTherminator::TTherminator(const TTherminator & therm) :
  59   TGenerator(therm),
  60   fCalka(0),
  61   fEvent(0),
  62   fPartDB(0)
  63 {
  64   // Copy constructor
  65   fPartDB = new ParticleDB();
  66 }
  67 TTherminator::~TTherminator()
  68 {
  69   // Destructor
  70   if (sRPInstance) delete sRPInstance;
  71   if (fEvent) delete fEvent;
  72   if (fCalka) delete fCalka;
  73   if (fPartDB) delete fPartDB;
  74 }
  75
  76 void TTherminator::ReadParameters()
  77 {
  78   // Read in global model parameters
  79   STR tModel;
  80   STR tTable;
  81
  82   try {
  83     sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
  84     tModel = sRPInstance->getPar("FreezeOutModel");
  85     if (tModel.Contains("SingleFreezeOut"))
  86       sModel = 0;
  87 /*MCH begin*/
  88     else if (tModel.Contains("Lhyquid3D"))
  89       sModel = 10;
  90     else if (tModel.Contains("Lhyquid2D"))
  91       sModel = 11;
  92 /*MCH end*/
  93     else if (tModel.Contains("BlastWaveVTDelay"))
  94       sModel = 6;
  95     else if (tModel.Contains("BlastWaveVT"))
  96       sModel = 2;
  97     else if (tModel.Contains("BlastWaveVLinearFormation"))
  98       sModel = 7;
  99     else if (tModel.Contains("BlastWaveVLinearDelay"))
 100       sModel = 8;
 101     else if (tModel.Contains("BlastWaveVLinear"))
 102       sModel = 4;
 103     else if (tModel.Contains("FiniteHydro"))
 104       sModel = 5;
 105     else {
 106       PRINT_MESSAGE("Unknown model type: " << tModel.Data());
 107       PRINT_MESSAGE("Please provide the proper name");
 108       exit(0);
 109     }
 110     if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
 111       sRunType = 4;
 112     else
 113       sRunType = 3;
 114     tTable = sRPInstance->getPar("TableType");
 115     if (tTable.Contains("Mathematica"))
 116       sTables = 0;
 117     else if (tTable.Contains("SHARE"))
 118       sTables = 1;
 119     else {
 120       PRINT_MESSAGE("Unknown table type: " << tTable.Data());
 121       exit(0);
 122     }
 123     sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
 124     fInputDir = sRPInstance->getPar("InputDirSHARE");
 125   }
 126   catch (STR tError) {
 127     PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
 128     PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
 129     exit(0);
 130   }
 131 }
 132
 133 void        TTherminator::Initialize(){
 134   // Initialize global variables
 135   // and particle and decay tables
 136   Parser       *tParser;
 137
 138   sRPFileName = "therminator.in";
 139   sRPInstance = new ReadPar(sRPFileName.Data());
 140
 141   ReadParameters();
 142
 143   tParser = new Parser(fPartDB);
 144
 145   tParser->ReadShare();
 146
 147   fCalka = new Integrator(sIntegrateSample);
 148   fCalka->ReadMultInteg(fPartDB);
 149
 150
 151   // Read in particle table from share
 152   // and add missing particle type to TDatabasePDG
 153
 154   char str[50];
 155   char str1[200];
 156
 157   //  ParticleType *tPartBuf;
 158
 159   TDatabasePDG *tInstance = TDatabasePDG::Instance();
 160   TParticlePDG *tParticleType;
 161
 162   //  AliWarning(Form("Reading particle types from particles.data"));
 163
 164   ifstream in((fInputDir+"/"+"particles.data").Data());
 165   //  ifstream in("particles.data");
 166
 167   int charge;
 168
 169   int number=0;
 170   if ((in) && (in.is_open()))
 171     {
 172       //START OF HEAD-LINE
 173       in.ignore(200,'\n');
 174       in.ignore(200,'\n');
 175       in.ignore(200,'\n');
 176       //END OF HEAD-LINE
 177
 178       while (in>>str)
 179         {
 180           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 181             {
 182               in.getline(str1,200);
 183               continue;
 184             }
 185           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
 186
 187           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
 188           number++;
 189           tParticleType = tInstance->GetParticle((int) mc);
 190           if (!tParticleType) {
 191             charge = 0;
 192             if (strstr(str, "plu")) charge = 1;
 193             if (strstr(str, "min")) charge = -1;
 194             if (strstr(str, "plb")) charge = -1;
 195             if (strstr(str, "mnb")) charge = 1;
 196             if (strstr(str, "plp")) charge = 2;
 197             if (strstr(str, "ppb")) charge = -2;
 198             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
 199             //      printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
 200             //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
 201             //      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));
 202
 203           }
 204         }
 205       in.close();
 206     }
 207 }
 208
 209 void        TTherminator::GenerateEvent()
 210 {
 211   // Run single event generation
 212   if ((sRunType == 3) || (sRunType == 4))
 213     {
 214       if (sRunType==4)
 215         fCalka->Randomize();
 216
 217       if (!fEvent) fEvent = new Event(fPartDB, fCalka);
 218       if (sRunType == 4) fEvent->Randomize();
 219
 220       fEvent->GenerateEvent();
 221       fEvent->DecayParticles();
 222       //      fEvent->WriteEvent(0);
 223     }
 224 }
 225
 226 Int_t       TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
 227 {
 228   // Import particles from a generated event into an external array
 229   const double kFmToGev = 0.197327;
 230
 231   if (particles == 0) return 0;
 232   TClonesArray &particlesR = *particles;
 233   particlesR.Clear();
 234   Int_t nump = 0;
 235   if (!fEvent) return 0;
 236   Int_t numpart = fEvent->GetParticleCount();
 237   printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
 238   for (Int_t iPart=0; iPart<numpart; iPart++) {
 239     Particle *tPart = fEvent->GetParticleOfCount(iPart);
 240     Int_t tFather;
 241     tFather = tPart->GetFather();
 242
 243     if (tFather> -1) {
 244       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
 245       mother->SetLastDaughter(iPart);
 246       if (mother->GetFirstDaughter()==-1)
 247         mother->SetFirstDaughter(iPart);
 248     }
 249     nump++;
 250     //    printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
 251     new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
 252                                       tFather, -1, -1, -1,
 253                                       tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
 254                                       tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
 255     particlesR[iPart]->SetUniqueID(iPart);
 256   }
 257
 258   return nump;
 259 }
 260
 261 TObjArray*  TTherminator::ImportParticles(Option_t */*option*/)
 262 {
 263   // Import generated particles into an internal array
 264   const double kFmToGev = 0.197327;
 265
 266   Int_t nump = 0;
 267   fParticles->Clear();
 268   if (!fEvent) return 0;
 269   Int_t numpart = fEvent->GetParticleCount();
 270   printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
 271   fEvent->InitParticleScan();
 272   for (Int_t iPart=0; iPart<numpart; iPart++) {
 273     Particle *tPart = fEvent->GetNextParticle();
 274     Int_t tFather;
 275     tFather = tPart->GetFather();
 276
 277     if (tFather> -1) {
 278       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
 279       mother->SetLastDaughter(iPart);
 280       if (mother->GetFirstDaughter()==-1)
 281         mother->SetFirstDaughter(iPart);
 282     }
 283     TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
 284                                  tFather, -1, -1, -1,
 285                                  tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
 286                                  tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
 287     p->SetUniqueID(iPart);
 288     fParticles->Add(p);
 289   }
 290   nump = fEvent->GetParticleCount();
 291
 292   return fParticles;
 293
 294 }