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   if (fCalka) delete fCalka;
  67   fCalka = new Integrator(*therm.fCalka);
  68   if (fEvent) delete fEvent;
  69   fEvent = new Event(*therm.fEvent);
  70   if (fPartDB) delete fPartDB;
  71   fPartDB = new ParticleDB();
  72 }
  73 TTherminator& TTherminator::operator=(const TTherminator & therm)
  74 {
  75   if (this != &therm) {
  76     fCalka = therm.fCalka;
  77     fEvent = therm.fEvent;
  78     delete fPartDB;
  79     fPartDB = new ParticleDB();
  80   }
  81
  82   return *this;
  83 }
  84
  85 TTherminator::~TTherminator()
  86 {
  87   // Destructor
  88   if (sRPInstance) delete sRPInstance;
  89   if (fEvent) delete fEvent;
  90   if (fCalka) delete fCalka;
  91   if (fPartDB) delete fPartDB;
  92 }
  93
  94 void TTherminator::ReadParameters()
  95 {
  96   // Read in global model parameters
  97   STR tModel;
  98   STR tTable;
  99
 100   try {
 101     sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
 102     tModel = sRPInstance->getPar("FreezeOutModel");
 103     if (tModel.Contains("SingleFreezeOut"))
 104       sModel = 0;
 105 /*MCH begin*/
 106     else if (tModel.Contains("Lhyquid3D"))
 107       sModel = 10;
 108     else if (tModel.Contains("Lhyquid2D"))
 109       sModel = 11;
 110 /*MCH end*/
 111     else if (tModel.Contains("BlastWaveVTDelay"))
 112       sModel = 6;
 113     else if (tModel.Contains("BlastWaveVT"))
 114       sModel = 2;
 115     else if (tModel.Contains("BlastWaveVLinearFormation"))
 116       sModel = 7;
 117     else if (tModel.Contains("BlastWaveVLinearDelay"))
 118       sModel = 8;
 119     else if (tModel.Contains("BlastWaveVLinear"))
 120       sModel = 4;
 121     else if (tModel.Contains("FiniteHydro"))
 122       sModel = 5;
 123     else {
 124       PRINT_MESSAGE("Unknown model type: " << tModel.Data());
 125       PRINT_MESSAGE("Please provide the proper name");
 126       exit(0);
 127     }
 128     if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
 129       sRunType = 4;
 130     else
 131       sRunType = 3;
 132     tTable = sRPInstance->getPar("TableType");
 133     if (tTable.Contains("Mathematica"))
 134       sTables = 0;
 135     else if (tTable.Contains("SHARE"))
 136       sTables = 1;
 137     else {
 138       PRINT_MESSAGE("Unknown table type: " << tTable.Data());
 139       exit(0);
 140     }
 141     sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
 142     fInputDir = sRPInstance->getPar("InputDirSHARE");
 143   }
 144   catch (STR tError) {
 145     PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
 146     PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
 147     exit(0);
 148   }
 149 }
 150
 151 void        TTherminator::Initialize(){
 152   // Initialize global variables
 153   // and particle and decay tables
 154   Parser       *tParser;
 155
 156   sRPFileName = "therminator.in";
 157   sRPInstance = new ReadPar(sRPFileName.Data());
 158
 159   ReadParameters();
 160
 161   tParser = new Parser(fPartDB);
 162
 163   tParser->ReadShare();
 164
 165   fCalka = new Integrator(sIntegrateSample);
 166   fCalka->ReadMultInteg(fPartDB);
 167
 168
 169   // Read in particle table from share
 170   // and add missing particle type to TDatabasePDG
 171
 172   char str[50];
 173   char str1[200];
 174
 175   //  ParticleType *tPartBuf;
 176
 177   TDatabasePDG *tInstance = TDatabasePDG::Instance();
 178   TParticlePDG *tParticleType;
 179
 180   //  AliWarning(Form("Reading particle types from particles.data"));
 181
 182   ifstream in((fInputDir+"/"+"particles.data").Data());
 183   //  ifstream in("particles.data");
 184
 185   int charge;
 186
 187   int number=0;
 188   if ((in) && (in.is_open()))
 189     {
 190       //START OF HEAD-LINE
 191       in.ignore(200,'\n');
 192       in.ignore(200,'\n');
 193       in.ignore(200,'\n');
 194       //END OF HEAD-LINE
 195
 196       while (in>>str)
 197         {
 198           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 199             {
 200               in.getline(str1,200);
 201               continue;
 202             }
 203           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
 204
 205           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
 206           number++;
 207           tParticleType = tInstance->GetParticle((int) mc);
 208           if (!tParticleType) {
 209             charge = 0;
 210             if (strstr(str, "plu")) charge = 1;
 211             if (strstr(str, "min")) charge = -1;
 212             if (strstr(str, "plb")) charge = -1;
 213             if (strstr(str, "mnb")) charge = 1;
 214             if (strstr(str, "plp")) charge = 2;
 215             if (strstr(str, "ppb")) charge = -2;
 216             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
 217             //      printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
 218             //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
 219             //      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));
 220
 221           }
 222         }
 223       in.close();
 224     }
 225
 226   delete tParser;
 227 }
 228
 229 void        TTherminator::GenerateEvent()
 230 {
 231   // Run single event generation
 232   if ((sRunType == 3) || (sRunType == 4))
 233     {
 234       if (sRunType==4)
 235         fCalka->Randomize();
 236
 237       if (!fEvent) fEvent = new Event(fPartDB, fCalka);
 238       if (sRunType == 4) fEvent->Randomize();
 239
 240       fEvent->GenerateEvent();
 241       fEvent->DecayParticles();
 242       //      fEvent->WriteEvent(0);
 243     }
 244 }
 245
 246 Int_t       TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
 247 {
 248   // Import particles from a generated event into an external array
 249   const double kFmToGev = 0.197327;
 250
 251   if (particles == 0) return 0;
 252   TClonesArray &particlesR = *particles;
 253   particlesR.Clear();
 254   Int_t nump = 0;
 255   if (!fEvent) return 0;
 256   Int_t numpart = fEvent->GetParticleCount();
 257   printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
 258   for (Int_t iPart=0; iPart<numpart; iPart++) {
 259     Particle *tPart = fEvent->GetParticleOfCount(iPart);
 260     Int_t tFather;
 261     tFather = tPart->GetFather();
 262
 263     if (tFather> -1) {
 264       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));
 265       mother->SetLastDaughter(iPart);
 266       if (mother->GetFirstDaughter()==-1)
 267         mother->SetFirstDaughter(iPart);
 268     }
 269     nump++;
 270     //    printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
 271     new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
 272                                       tFather, -1, -1, -1,
 273                                       tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
 274                                       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);
 275     particlesR[iPart]->SetUniqueID(iPart);
 276   }
 277
 278   return nump;
 279 }
 280
 281 TObjArray*  TTherminator::ImportParticles(Option_t */*option*/)
 282 {
 283   // Import generated particles into an internal array
 284   const double kFmToGev = 0.197327;
 285
 286   Int_t nump = 0;
 287   fParticles->Clear();
 288   if (!fEvent) return 0;
 289   Int_t numpart = fEvent->GetParticleCount();
 290   printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
 291   fEvent->InitParticleScan();
 292   for (Int_t iPart=0; iPart<numpart; iPart++) {
 293     Particle *tPart = fEvent->GetNextParticle();
 294     Int_t tFather;
 295     tFather = tPart->GetFather();
 296
 297     if (tFather> -1) {
 298       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));
 299       mother->SetLastDaughter(iPart);
 300       if (mother->GetFirstDaughter()==-1)
 301         mother->SetFirstDaughter(iPart);
 302     }
 303     TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
 304                                  tFather, -1, -1, -1,
 305                                  tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
 306                                  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);
 307     p->SetUniqueID(iPart);
 308     fParticles->Add(p);
 309   }
 310   nump = fEvent->GetParticleCount();
 311
 312   return fParticles;
 313
 314 }