TTherminator/Therminator/Event.cxx

   1 /******************************************************************************
   2  *                      T H E R M I N A T O R                                 *
   3  *                   THERMal heavy-IoN generATOR                              *
   4  *                           version 1.0                                      *
   5  *                                                                            *
   6  * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
   7  *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
   8  * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
   9  *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
  10  * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
  11  *                 Wojciech Florkowski                                        *
  12  *                                                                            *
  13  * For the detailed description of the program and furhter references         *
  14  * to the description of the model plesase refer to: nucl-th/0504047,         *
  15  * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
  16  *                                                                            *
  17  * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
  18  *                                                                            *
  19  * This code can be freely used and redistributed. However if you decide to   *
  20  * make modifications to the code, please contact the authors, especially     *
  21  * if you plan to publish the results obtained with such modified code.       *
  22  * Any publication of results obtained using this code must include the       *
  23  * reference to nucl-th/0504047 and the published version of it, when         *
  24  * available.                                                                 *
  25  *                                                                            *
  26  *****************************************************************************/
  27 #include <iostream>
  28 #include <TFile.h>
  29 #include <TH1D.h>
  30 #include "Particle.h"
  31 #include "Event.h"
  32 #include "Integrator.h"
  33 #include "ParticleDecayer.h"
  34 #include "DecayTable.h"
  35 #include "ReadPar.h"
  36
  37 extern ReadPar* sRPInstance;
  38
  39 Event::Event(ParticleDB *aDB, Integrator* aInteg) :
  40   mCurIter(),
  41   mScanStarted(0)
  42 {
  43   mDB = aDB;
  44   mInteg = aInteg;
  45   mRandom = new TRandom2();
  46 #ifdef _ROOT_4_
  47   mRandom->SetSeed2(31851, 14327);
  48 #else
  49   mRandom->SetSeed(31851);
  50 #endif
  51   mAverageMultiplicities.resize(mDB->GetParticleTypeCount());
  52   mMultiplicities.resize(mDB->GetParticleTypeCount());
  53   ReadParameters();
  54   ReadMultiplicities();
  55
  56
  57 }
  58
  59 Event::~Event()
  60 {
  61
  62   mParticles.clear();
  63   delete mRandom;
  64 }
  65
  66 void
  67 Event::GenerateEvent(int aSeed)
  68 {
  69   Particle **tPartBuf;
  70
  71   mParticles.clear();
  72 #ifdef _ROOT_4_
  73   if (aSeed) mRandom->SetSeed2(aSeed, (aSeed*2) % (7*11*23*31));
  74 #else
  75   if (aSeed) mRandom->SetSeed(aSeed);
  76 #endif
  77   GenerateMultiplicities();
  78   for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++) {
  79     mInteg->Generate(mDB->GetParticleType(tIter), mMultiplicities[tIter], &tPartBuf);
  80     for (int tBufIter=0; tBufIter<mMultiplicities[tIter]; tBufIter++) {
  81       AddParticle(tPartBuf[tBufIter]);
  82       delete tPartBuf[tBufIter];
  83     }
  84     free(tPartBuf);
  85   }
  86 }
  87
  88 int
  89 Event::GetParticleCount()
  90 {
  91   return mParticles.size();
  92 }
  93
  94 Particle*
  95 Event::GetNextParticle()
  96 {
  97   if (mScanStarted)
  98     InitParticleScan();
  99   if (mCurIter == mParticles.end())
 100     return 0;
 101
 102   Particle *tBuf;
 103   tBuf = &(*mCurIter);
 104   mCurIter++;
 105
 106   return tBuf;
 107 }
 108
 109 void
 110 Event::InitParticleScan()
 111 {
 112   mCurIter = mParticles.begin();
 113   mScanStarted = 1;
 114 }
 115
 116 Particle*
 117 Event::GetParticleOfCount(int aCount)
 118 {
 119   if (aCount == 0) {
 120     InitParticleScan();
 121   }
 122   else if (aCount == 1) {
 123     InitParticleScan();
 124     mCurIter++;
 125   }
 126   else {
 127     mCurIter--;
 128     mCurIter++;
 129     mCurIter++;
 130   }
 131
 132   if (mCurIter == mParticles.end())
 133     return 0;
 134
 135   Particle *tBuf;
 136   tBuf = &(*mCurIter);
 137   return tBuf;
 138 }
 139
 140 void
 141 Event::WriteEvent(int nEvent)
 142 {
 143   ParticleListIter tPLIter;
 144   tPLIter = mParticles.begin();
 145   int tCount = 0;
 146   ofstream os;
 147
 148   if (nEvent == 0) {
 149     try {
 150       os.open(sRPInstance->getPar("EventOutputFile").Data());
 151     }
 152     catch (STR e) {
 153       PRINT_MESSAGE("Very strange: No event output filename ???");
 154       exit(0);
 155     }
 156   }
 157   else {
 158     try {
 159       os.open(sRPInstance->getPar("EventOutputFile").Data(), ios::app);
 160     }
 161     catch (STR e) {
 162       PRINT_MESSAGE("Very strange: No event output filename ???");
 163       exit(0);
 164     }
 165   }
 166
 167
 168
 169   if (nEvent == 0)
 170     {
 171       os << "Therminator text output\nall_particles_p_x\nTherminator " << endl;
 172     }
 173   os << nEvent+1 << "   " << GetParticleCount() << "    " << "0.000" << "       " << "0.000" << endl;
 174
 175   while (tPLIter != mParticles.end())
 176     {
 177       os.flags(ios::right);
 178       os.precision(6);
 179       os.width(6);
 180       os << tCount << "   ";
 181       (*tPLIter).WriteParticle(&os);
 182       os << endl;
 183
 184       tCount++;
 185       tPLIter++;
 186     }
 187 }
 188
 189
 190 void
 191 Event::AddParticle(Particle* aParticle)
 192 {
 193   if (0) {
 194     ParticleListIter tPLIter;
 195     tPLIter = mParticles.begin();
 196     while (tPLIter != mParticles.end())
 197       {
 198         if (((*tPLIter).GetMass() == aParticle->GetMass()) &&
 199             ((*tPLIter).GetI3() == aParticle->GetI3()) &&
 200             ((*tPLIter).GetBarionN() == aParticle->GetBarionN()) &&
 201             ((*tPLIter).GetStrangeness() == aParticle->GetStrangeness())) {
 202           break; }
 203         if ((*tPLIter).GetMass() < aParticle->GetMass()) {
 204           break; }
 205         else if ((*tPLIter).GetMass() == aParticle->GetMass()) {
 206           if ((*tPLIter).GetI3() < aParticle->GetI3()) {
 207             break; }
 208           else if ((*tPLIter).GetI3() == aParticle->GetI3()) {
 209             if ((*tPLIter).GetBarionN() < aParticle->GetBarionN()) {
 210               break; }
 211             else if ((*tPLIter).GetBarionN() == aParticle->GetBarionN()) {
 212               if ((*tPLIter).GetStrangeness() < aParticle->GetStrangeness()) {
 213                 break; }
 214             }
 215           }
 216         }
 217         tPLIter++;
 218       }
 219     mParticles.insert(tPLIter, *aParticle);
 220   }
 221   else
 222     mParticles.push_back(*aParticle);
 223 }
 224
 225 void
 226 Event::ReadMultiplicities()
 227 {
 228   char   tName[200];
 229   double tMult;
 230
 231   char *tHash;
 232   int   multsize;
 233   char *tMultName = 0;
 234   ifstream *fin = NULL;
 235
 236   tHash = mInteg->ParameterHash();
 237   multsize = strlen(mFOHSlocation.Data()) + 25 + strlen(tHash);
 238   tMultName = (char *) malloc(sizeof(char) * multsize);
 239
 240   if (!tMultName) {
 241     printf("Cannot allocate memory!\n");
 242     exit(0);
 243   }
 244
 245   if (mFOHSlocation != "") {
 246     strcpy(tMultName, mFOHSlocation.Data());
 247     strcat(tMultName, "/");
 248     strcat(tMultName, "fmultiplicity_");
 249     strcat(tMultName, tHash);
 250     strcat(tMultName, ".txt");
 251     fin = new ifstream(tMultName);
 252   }
 253   else if (!((fin) && (fin->is_open()))) {
 254     strcpy(tMultName, "fmultiplicity_");
 255     strcat(tMultName, tHash);
 256     strcat(tMultName, ".txt");
 257     fin = new ifstream(tMultName);
 258   }
 259
 260 //   strcpy(tMultName, "fmultiplicity_");
 261 //   strcat(tMultName, tHash);
 262 //   strcat(tMultName, ".txt");
 263
 264 //  fin = new ifstream(tMultName);
 265   if ((fin) && (fin->is_open())) {
 266     PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
 267
 268     while (!fin->eof())
 269     {
 270       (*fin) >> tName >> tMult;
 271       PRINT_DEBUG_2(tName << " " <<  mDB->GetParticleTypeIndex(strdup(tName)) << " " << tMult);
 272       mAverageMultiplicities[mDB->GetParticleTypeIndex(strdup(tName))] = tMult;
 273     }
 274     fin->close();
 275   }
 276
 277   delete tMultName;
 278 }
 279
 280 void
 281 Event::GenerateMultiplicities()
 282 {
 283   for (unsigned int tIter=0; tIter<mAverageMultiplicities.size(); tIter++)
 284     mMultiplicities[tIter] = mRandom->Poisson(mAverageMultiplicities[tIter]);
 285 }
 286
 287 void
 288 Event::DecayParticles()
 289 {
 290   Particle *tPart1, *tPart2, *tPart3, *tFather;
 291   ParticleDecayer* tDecayer;
 292   int tCount = 0;
 293
 294   tDecayer = new ParticleDecayer(mDB);
 295   tDecayer->SeedSet(mRandom->Integer(100000000));
 296
 297   //  InitParticleScan();
 298   while ((tFather = GetParticleOfCount(tCount)))
 299     {
 300       if (tFather->GetParticleType()->GetGamma() >= 0.0) {
 301         if ((tFather->GetParticleType()->GetTable()) && (((DecayTable *) tFather->GetParticleType()->GetTable())->GetChannelCount()+1 > 0))
 302           {
 303             tDecayer->DecayParticle(tFather, &tPart1, &tPart2, &tPart3);
 304 #ifndef _RESCALE_CHANNELS_
 305             if (!tPart1) {
 306               tCount++;
 307               continue;
 308             }
 309 #endif
 310 #ifdef _NO_THREE_BODY_DECAYS_
 311             if (tPart3){
 312               tCount++;
 313               continue;
 314             }
 315 #endif
 316 #ifdef _OMIT_TWO_BODY_
 317             if (tPart1 && tPart2 && !tPart3){
 318               tCount++;
 319               continue;
 320             }
 321 #endif
 322
 323             tPart1->SetFather(tCount);
 324             tPart2->SetFather(tCount);
 325             if (tPart3)
 326               tPart3->SetFather(tCount);
 327
 328             AddParticle(tPart1);
 329             AddParticle(tPart2);
 330
 331             delete tPart1;
 332             delete tPart2;
 333
 334             if (tPart3) {
 335               AddParticle(tPart3);
 336               delete tPart3;
 337             }
 338           }
 339         else
 340           {
 341           }
 342       }
 343       tCount++;
 344     }
 345
 346   delete tDecayer;
 347 }
 348
 349 void
 350 Event::Randomize()
 351 {
 352   TDatime dat;
 353
 354 #ifdef _ROOT_4_
 355   mRandom->SetSeed2(dat.Get() / 2 * 3, dat.Get() / 11 * 9);
 356 #else
 357   mRandom->SetSeed(dat.Get() / 2 * 3);
 358 #endif
 359 }
 360
 361 void
 362 Event::ReadParameters()
 363 {
 364   STR tEvFile;
 365   STR tUseNegativeBinomial;
 366
 367   try {
 368     tEvFile = sRPInstance->getPar("EventOutputFile");
 369   }
 370   catch (STR e) {
 371     PRINT_DEBUG_1("Event::ReadParameters - Caught exception " << e);
 372     PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
 373     exit(0);
 374   }
 375   try {
 376     tUseNegativeBinomial = sRPInstance->getPar("MultiplicityDistribution");
 377     if (tUseNegativeBinomial.Contains("NegativeBinomial"))
 378       mNegBin = 1;
 379     else
 380       mNegBin = 0;
 381   }
 382   catch (STR e) {
 383     PRINT_MESSAGE("Using default multiplicty distribution: Poissonian");
 384     mNegBin = 0;
 385   }
 386   try {
 387     mFOHSlocation = sRPInstance->getPar("FOHSLocation");
 388   }
 389   catch (STR e) {
 390   }
 391
 392 }
 393