TUHKMgen/UHKM/RunHadronSource.cxx

   1 /*******************************************************************************
   2  *                                                                             *
   3  *    HYDJET++ , event generator under the ROOT FRAMEWORK for simulation of    *
   4  *    relativistic heavy ion AA collisions as the superposition of soft,       *
   5  *    hydro-type state and hard, multi-parton state.                           *
   6  *                                                                             *
   7  *     The main routine is written in the object-oriented C++ language         *
   8  *     under the ROOT environment. The hard, multi-partonic part of            *
   9  *     HYDJET++ event is identical to the hard part of Fortran-written         *
  10  *     HYDJET (PYTHIA6.4xx + PYQUEN1.5) and is included in the generator       *
  11  *     structure as the separate directory. The soft part of HYDJET++          *
  12  *     event represents the "thermal" hadronic state obtained with the         *
  13  *     parameterization Bjorken-like of freeze-out hypersurface and            *
  14  *     includes longitudinal, radial and elliptic flow effects and             *
  15  *     decays of hadronic resonances. The corresponding fast                   *
  16  *     Monte-Carlo simulation procedure (C++ code) FAST MC is adapted.         *
  17  * --------------------------------------------------------------              *
  18  *     Web-page:                                                               *
  19  *    http://cern.ch/lokhtin/hydjet++                                          *
  20  *     --------------------------------------------------------------          *
  21  *                                                                             *                                                                             *
  22  *                                                                             *
  23  * This program is a free software; you can use and redistribute it freely.    *
  24  * Any publication of results obtained using this code must reference          *
  25  *                                                                             *
  26  *                                                                             *
  27  *                                                                             *
  28  *      Main reference for HYDJET++:                                           *
  29  *     I.P. Lokhtin, L.V. Malinina, S.V. Petrushanko, A.M. Snigirev,           *
  30  *     I. Arsene, K. Tywoniuk, submitted to Comp. Phys. Comm.                  *
  31  *                                                                             *
  32  *     Reference for HYDJET and PYQUEN:                                        *
  33  *     I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211;             *
  34  *     http://cern.ch/lokhtin/hydro/hydjet.html                                *
  35  *     http://cern.ch/lokhtin/pyquen.                                          *
  36  *                                                                             *
  37  *     Reference for PYTHIA6.4:                                                *
  38  *     T.Sjostrand, S. Mrenna and P. Skands, JHEP05 (2006) 026;                *
  39  *     http://home.thep.lu.se/~torbjorn/Pythia.html.                           *
  40  *                                                                             *
  41  *     References for FAST MC:                                                 *
  42  *     N.S. Amelin, R. Lednicky, T.A. Pocheptsov, I.P. Lokhtin,                *
  43  *     L.V. Malinina, A.M. Snigirev, Iu.A. Karpenko and Yu.M. Sinyukov,        *
  44  *     Phys. Rev. C 74 (2006) 064901;                                          *
  45  *     N.S. Amelin, I. Arsene, L. Bravina, Iu.A. Karpenko, R. Lednicky,        *
  46  *     I.P. Lokhtin, L.V. Malinina, A.M. Snigirev and Yu.M. Sinyukov,          *
  47  *     Phys. Rev. C 77 (2008) 014903;                                          *
  48  *     http://uhkm.jinr.ru.                                                    *
  49  *                                                                             *
  50  *     Reference for nuclear shadowing model:                                  *
  51  *     K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and                   *
  52  *     E. Zabrodin, Phys. Lett. B 657 (2007) 170.                              *
  53  *                                                                             *
  54  *       version 2.0:                                                          *
  55  *                                                                             *
  56  *     Igor Lokhtin, SINP MSU, Moscow, RU                                      *
  57  *     e-mail: Igor.Lokhtin@cern.ch                                            *
  58  *                                                                             *
  59  *     Ludmila Malinina, SINP MSU, Moscow, RU                                  *
  60  *     e-mail: malinina@lav01.sinp.msu.ru                                      *
  61  *                                                                             *
  62  *******************************************************************************/
  63
  64 #include <iostream>
  65 #include <fstream>
  66 #include <vector>
  67 #include <time.h>
  68
  69 #include <TNtuple.h>
  70 #include <TError.h>
  71 #include <TTree.h>
  72 #include <TFile.h>
  73
  74 #include "InitialState.h"
  75 #include "InitialStateHydjet.h"
  76
  77 #include <TRandom.h>
  78 #include "Particle.h"
  79
  80
  81 #include "HYJET_COMMONS.h"
  82 extern HYIPARCommon HYIPAR;
  83 extern HYFPARCommon HYFPAR;
  84 extern HYJPARCommon HYJPAR;
  85 extern HYPARTCommon HYPART;
  86 extern SERVICECommon SERVICE;
  87
  88
  89 //Main program:
  90 //reads input parameters from file "RunInputHydjet" ;
  91 //calculates particle densities and average initial multiplicities and writes them
  92 //in output file "multiplicities.txt";
  93 //creates trees (tree with direct hadrons and hadrons after resonance decays)
  94 //with space-time and momentum-energy information of produced hadrons;
  95 //writes trees in file "RunOutput.root".
  96
  97
  98 int main()
  99
 100 {
 101
 102   clock_t start;
 103   start = clock();
 104
 105 //new
 106   time_t  now;
 107   struct tm  *ts;
 108   char       buf[80];
 109
 110  // Get the current time
 111    time(&now);
 112
 113  // Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
 114     ts = localtime(&now);
 115     strftime(buf, sizeof(buf), "%a %Y-%m-%d %H:%M:%S %Z", ts);
 116     printf("%s\n", buf);
 117
 118
 119
 120
 121   TFile *outputFile=new TFile("RunOutput.root", "RECREATE");
 122
 123   //SET MAXIMAl VALUE OF PARTICLE MULTIPLICITY!!!
 124   const Int_t kMax = 500000;
 125
 126   Int_t npart;
 127   //define event number
 128   Int_t nev;
 129
 130 //define event characteristics:
 131 //total event multiplicity, number of produced hadrons in hard part/soft part
 132   Int_t Ntot, Npyt, Nhyd;
 133
 134 // number of jets, number of binary collisions, number of participants :
 135   Int_t Njet, Nbcol, Npart;
 136
 137  //impact parameter
 138   Float_t Bgen, Sigin, Sigjet;
 139
 140   //define hadron characteristic vectors
 141   std::vector<Int_t> pdg(kMax); //pdg encodings
 142   std::vector<Int_t> Mpdg(kMax);//pdg encodings for mother hadrons
 143   std::vector<Int_t> type(kMax);//type of particle: 0-from hydro or decays, 1111 -from jets
 144   std::vector<Float_t> Px(kMax);//x-hadron momentum component,[GeV/c]
 145   std::vector<Float_t> Py(kMax);//y-hadron momentum component,[GeV/c]
 146   std::vector<Float_t> Pz(kMax);//z-hadron momentum component,[GeV/c]
 147   std::vector<Float_t> E(kMax); //hadron total energy,[GeV]
 148   std::vector<Float_t> X(kMax);//x-hadron coordinate component,[fm]
 149   std::vector<Float_t> Y(kMax);//y-hadron coordinate component,[fm]
 150   std::vector<Float_t> Z(kMax);//z-hadron coordinate component,[fm]
 151   std::vector<Float_t> T(kMax);//hadron time,[fm/c]
 152
 153
 154
 155   TTree *ti=new TTree("ti","Initial");
 156
 157   ti->Branch("nev",&nev,"nev/I");
 158   ti->Branch("Bgen",&Bgen,"Bgen/F");
 159   ti->Branch("Sigin",&Sigin,"Sigin/F");
 160   ti->Branch("Sigjet",&Sigjet,"Sigjet/F");
 161   ti->Branch("Ntot",&Ntot,"Ntot/I");
 162   ti->Branch("Nhyd",&Nhyd,"Nhyd/I");
 163   ti->Branch("Npyt",&Npyt,"Npyt/I");
 164   ti->Branch("Njet",&Njet,"Njet/I");
 165   ti->Branch("Nbcol",&Nbcol,"Nbcol/I");
 166   ti->Branch("Npart",&Npart,"Npart/I");
 167   ti->Branch("Px",&Px[0],"Px[Ntot]/F");
 168   ti->Branch("Py",&Py[0],"Py[Ntot]/F");
 169   ti->Branch("Pz",&Pz[0],"Pz[Ntot]/F");
 170   ti->Branch("E",&E[0],"E[Ntot]/F");
 171   ti->Branch("X",&X[0],"X[Ntot]/F");
 172   ti->Branch("Y",&Y[0],"Y[Ntot]/F");
 173   ti->Branch("Z",&Z[0],"Z[Ntot]/F");
 174   ti->Branch("T",&T[0],"T[Ntot]/F");
 175   ti->Branch("pdg",&pdg[0],"pdg[Ntot]/I");
 176   ti->Branch("Mpdg",&Mpdg[0],"Mpdg[Ntot]/I");
 177
 178   TTree *td=new TTree("td","After decays");
 179   td->Branch("nev",&nev,"nev/I");
 180   td->Branch("Bgen",&Bgen,"Bgen/F");
 181   td->Branch("Sigin",&Sigin,"Sigin/F");
 182   td->Branch("Sigjet",&Sigjet,"Sigjet/F");
 183   td->Branch("Ntot",&Ntot,"Ntot/I");
 184   td->Branch("Nhyd",&Nhyd,"Nhyd/I");
 185   td->Branch("Npyt",&Npyt,"Npyt/I");
 186   td->Branch("Njet",&Njet,"Njet/I");
 187   td->Branch("Nbcol",&Nbcol,"Nbcol/I");
 188   td->Branch("Npart",&Npart,"Npart/I");
 189   td->Branch("Px",&Px[0],"Px[Ntot]/F");
 190   td->Branch("Py",&Py[0],"Py[Ntot]/F");
 191   td->Branch("Pz",&Pz[0],"Pz[Ntot]/F");
 192   td->Branch("E",&E[0],"E[Ntot]/F");
 193   td->Branch("X",&X[0],"X[Ntot]/F");
 194   td->Branch("Y",&Y[0],"Y[Ntot]/F");
 195   td->Branch("Z",&Z[0],"Z[Ntot]/F");
 196   td->Branch("T",&T[0],"T[Ntot]/F");
 197   td->Branch("pdg",&pdg[0],"pdg[Ntot]/I");
 198   td->Branch("Mpdg",&Mpdg[0],"Mpdg[Ntot]/I");
 199   td->Branch("type",&type[0],"type[Ntot]/I");
 200
 201
 202   InitialState *FASTMC;
 203     FASTMC = new InitialStateHydjet();
 204
 205
 206   if(!FASTMC->ReadParams()) {
 207     Error("RunHadronSource::main", "No initial model parameters found!!\n");
 208     return 0;
 209   }
 210
 211   if(!FASTMC->MultIni()) {
 212     Error("RunHadronSource::main", "Initial multiplicities are zero!!\n");
 213     return 0;
 214   }
 215
 216   ParticleAllocator allocator;
 217   List_t source;
 218   List_t secondaries;
 219   std::cout << "Generating " << FASTMC->GetNev() << " events" << std::endl;
 220   std::cout << "Starting the event loop" << std::endl;
 221
 222   // Set Random Number seed
 223   Int_t sseed =0;               //Set 0 to use the current time
 224   gRandom->SetSeed(sseed);
 225   std::cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<std::endl;
 226
 227   // Loop over events
 228   for(Int_t ev = 0; ev < FASTMC->GetNev(); ++ev) {
 229     nev = ev;
 230     // Initialize the source
 231     FASTMC->Initialize(source, allocator);
 232
 233       Npart = HYFPAR.npart;
 234       Bgen = HYFPAR.bgen;
 235       Njet = HYJPAR.njet;
 236       Nbcol = HYFPAR.nbcol;
 237       if(ev==0) {
 238             Sigin=HYJPAR.sigin;
 239             Sigjet=HYJPAR.sigjet;
 240          }
 241      std::cout<<"in RunHadronSource: ev"<<ev<<" Njet "<<Njet<<" Nbcol "<<Nbcol<<" Npart "<<Npart<<std::endl;
 242
 243     if(source.empty()) {
 244       Error("RunHadronSource::main", "Source is not initialized!!");
 245       return 0;
 246     }
 247
 248     // Run the decays
 249     if(FASTMC->GetTime() >= 0.)
 250       FASTMC->Evolve(source, secondaries, allocator, FASTMC->GetWeakDecayLimit());
 251
 252     std::cout << "event #" << ev << "\r" << std::flush;
 253     Ntot = 0;
 254     Npyt = 0;
 255     Nhyd = 0;
 256     Ntot = 0;
 257
 258     LPIT_t it;
 259     LPIT_t e;
 260
 261     // Fill the source tree
 262
 263     Ntot = 0; Nhyd=0; Npyt=0;
 264     for(it = source.begin(), e = source.end(); it != e; ++it) {
 265       TVector3 pos(it->Pos().Vect());
 266       TVector3 mom(it->Mom().Vect());
 267       Float_t m1 = it->TableMass();
 268       pdg[Ntot] = it->Encoding();
 269       Mpdg[Ntot] = -1;
 270       Px[Ntot] = mom[0];
 271       Py[Ntot] = mom[1];
 272       Pz[Ntot] = mom[2];
 273       E[Ntot] = TMath::Sqrt(mom.Mag2() + m1*m1);
 274       X[Ntot] = pos[0];
 275       Y[Ntot] = pos[1];
 276       Z[Ntot] = pos[2];
 277       T[Ntot] = it->T();
 278       type[Ntot] = it->GetType();
 279       if(type[Ntot]==0)Nhyd++;
 280       if(type[Ntot]==1)Npyt++;
 281       Ntot++;
 282       if(Ntot > kMax)
 283         Error("in main:", "Ntot is too large %d", Ntot);
 284     }
 285     ti->Fill();
 286
 287     std::cout<<"ti Ntot= " <<Ntot<<" Npyt= "<<Npyt<<" Nhyd=  "<<Nhyd<<std::endl;
 288
 289     // Fill the decayed tree
 290     Ntot = 0; Nhyd=0; Npyt=0;
 291     for(it = secondaries.begin(), e = secondaries.end(); it != e; ++it) {
 292       TVector3 pos(it->Pos().Vect());
 293       TVector3 mom(it->Mom().Vect());
 294       Float_t m1 = it->TableMass();
 295       pdg[Ntot] = it->Encoding();
 296       Mpdg[Ntot] = it->GetLastMotherPdg();
 297       Px[Ntot] = mom[0];
 298       Py[Ntot] = mom[1];
 299       Pz[Ntot] = mom[2];
 300       E[Ntot] =  TMath::Sqrt(mom.Mag2() + m1*m1);
 301       X[Ntot] = pos[0];
 302       Y[Ntot] = pos[1];
 303       Z[Ntot] = pos[2];
 304       T[Ntot] = it->T();
 305       type[Ntot] = it->GetType();
 306       if(type[Ntot]==0)Nhyd++;
 307       if(type[Ntot]==1)Npyt++;
 308       Ntot++;
 309       if(Ntot > kMax)
 310         Error("in main:", "Ntot is too large %d", Ntot);
 311     }
 312     td->Fill();
 313
 314     std::cout<<"td Ntot= " <<Ntot<<" Npyt= "<<Npyt<<" Nhyd=  "<<Nhyd<<std::endl;
 315
 316     allocator.FreeList(source);
 317     allocator.FreeList(secondaries);
 318   }
 319
 320
 321
 322
 323   // Close the output file by getting it from the tree object to avoid crashes when
 324   // the output file is automatically switched by root.
 325   TFile *saveFile = ti->GetCurrentFile();
 326   saveFile->cd();
 327   ti->Write();
 328   td->Write();
 329   saveFile->Close();
 330   //  ti->GetCurrentFile()->Close();
 331
 332   clock_t stop;
 333   stop = clock();
 334   std::cout << "*********************************************" << std::endl;
 335   std::cout << "Execution time: " << (stop - start)/CLOCKS_PER_SEC << " seconds" << std::endl;
 336   std::cout << "*********************************************" << std::endl;
 337
 338
 339 //new
 340   time_t  now1;
 341   struct tm  *ts1;
 342   char       buf1[80];
 343
 344  // Get the current time
 345    time(&now1);
 346
 347  // Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
 348     ts1 = localtime(&now1);
 349     strftime(buf1, sizeof(buf1), "%a %Y-%m-%d %H:%M:%S %Z", ts1);
 350     printf("%s\n", buf1);
 351
 352
 353
 354
 355
 356   //return 0;
 357 }